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1.  Objectives 


This  program  was  intended  to  develop  theory,  a  broad  range  of  algorithms  and  proof  of  concept  in  detection, 
estimation,  and  reconstruction  of  objects  embedded  in  turbid  media,  which  hamper  visibility  (such  as  fog, 
clouds,  smoke,  and  other  types  of  aerosol  particles).  In  particular,  the  program  addresses  the  following  eight 
specific  areas: 


1.  We  exploit  recent  advances  in  the  physical  design  of  fast  optical  systems  which  enable  active  imaging 
and  ranging  with  ’’ballistic”  light.  In  this  modality,  fast  bursts  of  optical  energy  are  transmitted  into  a 
medium,  and  the  ballistic  component  of  light  (which  travels  through  a  medium  with  minimal  diffusive 
distortion)  is  detected  in  either  backscatter  from  the  target,  or  in  trans-illuminated  form  through  the 
medium. 

2.  We  simultaneously  optimize  the  shape  and  duration  of  the  optical  pulses  and  the  design  of  the  optical 
detectors  to  achieve  maximum  signal  detectability.  The  currently  existing  ballistic  imaging  practice  has 
been  based  largely  on  ad-hoc  choices  of  the  optical  waveforms,  and  static  detector  designs  based  on 
time-gating. 

3.  We  study  the  ballistic  imaging  problem  as  (1)  a  detection  problem  where  the  simple  presence  or 
absence  of  some  target  and  its  rough  dimensions  are  sought;  and  (2)  a  reconstruction  problem,  where 
a  full  image  of  the  object  is  to  be  obtained  when  possible.  The  first  approach  is  based  on  statistical 
detection  theory,  whereas  the  second  is  a  statistical  estimation  problem. 
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4.  We  study  two  complementary  approaches  for  addressing  the  detection  problem  and  improving  the 
baseline  performance.  First,  we  formulate  and  use  a  forward  model  of  the  ballistic  imaging  process  to 
enable  the  determination  of  the  optimum  signaling  strategies.  Second,  to  improve  the  detection 
performance,  we  adopt  an  active  heterodyne  signaling/sensing  which  uses  temporal  modulation  of  light 
to  significantly  improve  the  detector  sensitivity  by  mitigating  the  ill  effects  of  amplifier  noise. 

5.  We  exploit  temporal,  spatial,  wavelength  diversity  and  coding  for  ballistic  imaging.  Most,  if  not  all, 
ballistic  imaging  methods  described  in  the  literature  employ  a  raster-scanned  point-by-point  imaging 
system.  Such  approaches  are  inherently  limited  as  they  do  not  take  full  advantage  of  the  degrees  of 
freedom  in  the  signal  space.  We  study  the  array  imaging  framework  where  multiple  emitters  will 
transmit  coordinated  pulses  of  light;  and  for  their  part,  a  collection  of  photon-detecting  elements  will 
gather  the  received  data  and  compute  a  resultant  image.  In  addition  to  this  spatial  component,  and  the 
temporal  modulations  introduced  for  the  heterodyne  signaling  framework,  the  wavelength  of  light  may 
also  be  used  to  advantage  and  adapted.  Furthermore,  there  is  the  possibility  of  analyzing  various  ways 
of  coding  these  pulses  (e.g.  transmitting  sequences  of  pulses  as  is  done  in  Ultra-Wide-Band 
communications). 

6.  We  study  novel  multi-scale  approaches  to  the  photon-limited  image  reconstruction  problem.  Recently 
developed  methodology  in  this  area  will  enable  imaging  through  turbid  media  at  very  low  average 
photon  counts.  Meanwhile,  adaptive  sampling  strategies,  motivated  by  the  same  multi-scale  models, 
have  also  been  considered.  These  will  allow  for  a  coarse-to-fine  sampling  strategy  which  will  bootstrap 
few  coarse  (but  less  noisy)  samples,  to  guide  the  selection  of  higher  resolution  samples  to  form  a 
complete  image. 

7.  Theoretical  models  for  the  detection  and  estimation  problems  developed  above  are  subjected  to 
sensitivity  analysis,  hence  yielding  fundamental  performance  bounds  in  the  form  of  (1)  analytical 
relationships  between  minimum  required  SNR’s  to  achieve  a  particular  detection  power,  and  (2)  mean- 
squared  error  bounds  that  will  demonstrate  how  the  performance  of  the  image  estimation  problem  will 
depend  upon  the  underlying  physical  and  sensor  parameters.  Relationships  thus  derived  can  then  be 
used  to  optimize  system  performance  (e.g.  by  minimizing  the  Cramer-Rao  bound  for  instance  for  an 
estimation  problem.) 

8.  Proof  of  concept  is  carried  out  in  both  the  desk-top  and  the  bench-top.  We  build  simulated  models  of 
the  forward  ballistic  light  transport  process  and  test  the  signaling  strategies  and  detection  algorithms 
with  this  model. 


2 


2.  Status  of  Effort 


The  work  during  the  past  15  month  has  culminated  in  several  peer  reviewed  publications,  which  are  attached  to 
this  report.  The  topics  we  have  addressed  include  the  following: 

1 .  Collection  of  existing  data  for  light  propagation  through  turbid  media  which  enabled  us  to  move  forward 
with  the  analytical  work  without  having  to  wait  for  the  experimental  setup  to  be  completed.  As  will  be 
explained  later  in  this  document,  aside  from  the  data  generated  at  our  UCSC  lab,  we  have 
experimented  on  several  data  sets  generated  by  the  leading  physics  research  groups  in  the  area  of 
ultrafast  ballistic  imaging  (City  University  of  New  York  and  Universita  della  Tuscia,  Italy). 

2.  Forward  Modeling  and  Simulation  of  ballistic  and  diffuse  data.  A  forward  model  of  light  propagation 
through  turbid  media  allows  us  to  formulate  analytical  detection  and  estimation  problems  and  study  the 
performance  bounds.  We  have  combed  the  literature,  and  identified  the  most  appropriate  and 
parsimonious  models  which  can  be  used  to  describe  light  propagation  through  turbid  media.  We  have 
validated  these  models  using  statistical  means  (Adjusted  R2  goodness  of  fit)  against  real  data. 

3.  UCSC  laboratory  setup  and  data  collection  are  completed  which  enabled  us  to  produce  data  under 
controlled  conditions  to  aid  model  development  and  verification  of  the  performance  of  proposed 
analytical  algorithms.  UCSC  laboratory  is  developed  in  order  to  be  flexible  for  active  waveform  design. 
Single  point  optical  autocorrelation,  and  repeatable  Ballistic  Imaging  experiments  have  already  been 
demonstrated  in  this  laboratory. 

4.  Development  of  detection  algorithms  based  upon  ballistic  and  diffuse  data  is  completed.  In  trans¬ 
illumination  mode,  such  algorithms  enable  the  detection  of  objects  present  between  a  transmitting  light 
source  and  a  detector.  The  output  of  this  process  is  a  “binary”  image,  indicating  the  presence  or 
absence  of  obstructing  objects.  For  instance,  looking  down  a  foggy  road,  we  can  anticipate  the 
presence  of  enemy  soldiers  or  vehicles.  We  have  derived  ROC  curves  for  detection  of  a  canonical 
object  using  current  state  of  the  art  (time-gating)  methods.  We  showed  improvement  upon  these  ROC 
curves  by  developing  appropriate  optimal  statistical  signal  processing  methods. 

5.  Development  of  image  reconstruction  and  denoising  algorithms  based  upon  ballistic  and  diffuse  data  is 
completed.  Image  reconstruction  allows  for  computation  of  a  more  sophisticated  output.  Namely  the 
output  of  this  stage  is  a  gray-level  image  (or  a  sequence  of  images)  that  indicate  not  only  the  presence 
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or  absence  of  objects,  but  also  the  relative  reflectivity  (in  reflection  mode)  or  opaqueness  (in  trans¬ 
illumination  mode)  of  objects  of  interest  in  the  turbid  medium.  Using  the  measure  of  relative  PSNR 
(Peak  Signal  to  Noise  Ratio),  compared  to  the  state  of  the  art  (time-gating)  image  capture,  images 
reconstructed  using  methods  developed  by  Milanfar  et  al.  and  Nowak  et  al.  showed  improvement  of 
about  10  dB  relative  to  direct  ballistic  image  reconstruction  procedures  previously  proposed  in  the 
literature. 

6.  Development  of  performance  characteristics/bounds  on  above  algorithms  based  upon  the  forward 

model  is  completed.  Performance  analysis  and  bounds  enabled  us  to  understand  the  limitations 
inherent  in  the  solution  of  the  detection/estimation  problems.  Understanding  these  limitations  will  help 
in  defining  what  is  possible  and  impossible.  We  have  already  provided  analytical  formulae  for 

predicting  bias,  variance  (or  lower  bounds  to  these)  for  detection/estimation  algorithms. 

7.  Development  of  waveform  design  methodologies  to  maximize  the  performance  of  the  above 
detection/reconstruction  algorithms.  Using  the  performance  metrics  defined  above,  better  design  of 
waveforms  and  corresponding  light  detectors  will  enable  improve  visibility  in  turbid  media.  Some 
preliminary  theoretical  work  has  been  already  conducted  and  further  work  is  underway.. 


3.  Accomplishments 

Below  we  present  the  abstracts  of  our  key  recent  accomplishments  which  detail  the  progress  made  to  date 
on  problems  of  concern  to  this  program. 

•  We  have  investigated  and  validated  a  computationally  efficient,  yet  sufficiently  accurate  mathematical 
model  for  light  propagation  in  homogeneous  (target  absent)  and  inhomogeneous  (target  present)  turbid 
media. 

•  We  have  derived  the  fundamental  limits  on  the  accuracy  of  the  estimated  parameters  of  the  said 
mathematical  model  of  an  unknown  turbid  medium.  This  study  also  guides  us  toward  the  most  efficient 
experiments  (with  respect  to  both  time/accuracy)  for  calibrating  the  model  parameters  of  the  unknown 
turbid  medium  as  well  as  the  optical  properties  of  the  target  (which  can  be  used  to  identify/categorize 
it).  Our  simulation  results  show  that  for  the  medium  of  most  interest,  namely  heavy  fog;  optical 
parameters  can  be  estimated  with  very  high  accuracy.  These  experiments  are  equally  valid  in 
assessing  the  effectiveness  of  our  imaging  techniques  in  virtually  any  unknown  turbid  medium. 
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•  We  used  the  said  model  to  derive  optimal  statistical  tests  for  detecting  objects  hidden  in  turbid  media. 
Performance  analysis  was  carried  out  by  computing  ROC  curves  for  the  proposed  optimal  tests, 
showing  that  by  considering  only  the  ballistic  photons,  we  are  able  to  detect  opaque  objects  hidden  in 
heavy  fog  in  the  range  of  approximately  380  meters  (i.e.  30  mean-free  paths  or  MFP's).  Detection  rate 
of  the  semi  transparent  objects  are  shown  to  be  slightly  less  than  this  distance. 

•  To  further  improve  the  detection  rate  of  the  aforementioned  single  pixel  optimal  detectors  and 
penetrate  longer  distances  in  turbid  media,  we  exploited  sampling  at  a  diversity  of  locations  in  space. 
We  then  developed  an  algorithm  based  upon  the  Generalized  Likelihood  Ratio  Test  (GLRT)  framework 
which  takes  advantage  of  the  spatial  correlation  of  nearby  samples.  In  several  experiments,  we 
showed  that  objects  of  different  size  and  shape  that  are  completely  unrecognizable  using  the  common 
single  pixel  detection  techniques,  are  detectable  with  very  high  accuracy  using  the  said  multi-scale 
technique. 

•  Considering  the  effect  of  diffraction  in  turbid  media,  we  have  derived  a  lower  bound  on  the  achievable 
spatial  resolution  of  the  proposed  imaging  systems.  A  basic  example  of  such  study  is  that  it  is 
theoretically  possible  to  resolve,  with  high  accuracy,  opaque  objects  of  radius  on  the  order  of  6mm  in 
heavy  fog  at  a  distance  of  200  meters. 

•  To  penetrate  farther  in  turbid  media,  and  more  importantly  to  reduce  the  image  scanning  time,  we  have 
studied  and  experimented  on  several  adaptive  sampling  techniques.  The  experiments  preformed  using 
our  ballistic  imaging  laser  scanner  at  UCSC,  demonstrate  that  we  may  reduce  the  scanning  (imaging) 
time  by  as  much  as  a  factor  of  60  (depending  on  the  shape/size  of  the  target),  which  is  extremely 
important  for  imaging  in  hostile  battlefield  environments.  A  theoretic  study  on  the  expected  gain  (in 
Mean  Square  Error  term)  of  our  proposed  adaptive  sampling  technique  is  also  derived  and  presented. 

•  Since  the  output  of  adaptive  scans  are  in  general  irregularly  sampled  images,  we  have  developed 
optimal  image  reconstruction  algorithms  to  effectively  reconstruct  a  high  resolution  image  from  a  set  of 
irregularly  sampled,  blurry,  noisy  images.  Moreover,  our  image  reconstruction  technique  is  especially 
useful  as  the  proposed  time  resolved  holistic  imaging  system  exploits  diffused  photons  as  well  as  the 
ballistic  ones. 

•  We  established  the  performance  tradeoffs  between  transillumination  (TRT)  imaging  and  conventional 
(un-gated)  imaging.  On  the  one  hand,  time-gated,  First  arrival  photons  are  unscattered  and  therefore 
provide  very  high  spatial  resolution.  But,  very  few  photons  arrive  at  the  detector  without  scattering, 
effectively  resulting  in  a  very  low  SNR.  On  the  other  hand,  conventional  (un-gated)  imaging  is  based  on 
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all  photons  (scattered  and  unscattered),  resulting  in  lower  spatial  resolution,  but  higher  SNR  (due  to 
the  large  number  of  photons).  This  tradeoff  was  analyzed  using  a  decision-theoretic  approach  to 
ascertain  bounds  on  the  minimum  resolvable  occluding  object  size  with  and  without  time-gated  photon 
acquisition.  The  theoretical  predictions  are  validated  through  a  realistic  simulation  of  tumors  in  breast 
tissue.  Our  mathematical  analysis  clearly  indicates  in  which  regimes  (of  turbidity  and  scattering)  TRT 
imaging  outperforms  conventional  imaging.  These  theoretical  predictions  match  very  well  with 
experimental  results.  We  also  developed  a  novel  algorithm  for  TRT  image  reconstruction,  based  on 
the  principle  of  Maximum  Likelihood  Estimation.  The  innovative  algorith  combines  multiple  snapshot 
observations  of  time-gated  photons  to  reconstruct  a  high-resolution,  low-noise  image.  In  effect,  the 
algorithm  allows  us  to  obtain  the  "best-of-both-worlds"  by  combining  the  spatial  resolution  of  the  first- 
arrival  photons,  with  the  higher  SNR  provided  by  scattered  photons. 

•  We  set  up  a  ballistic  imaging  laboratory  at  UCSC,  enabling  us  to  capture  ballistic  images  with  what  can 
be  thought  of  as  the  "optimal"  (500fs)  time  resolution,  as  this  is  virtually  the  resolution  of  the  generated 
ballistic  pulses. 


4.  Personnel  Supported 

Peyman  Milanfar,  PI 

Benjamin  Friedlander,  Co-Pi 

Ali  Shakouri,  Co-Pi 

Michael  Isaacson,  Co-Pi 

Robert  Nowak,  (Univ.  of  Wisconsin)  Co-Pi 

Jonathan  Saint  Clair,  (Boeing  Co.)  Co-Pi 

Post-docs: 

Sina  Farsiu 
James  Christofferson 

Graduate  Students: 

Brian  Eriksson 
Priyam  Chatterjee 
Hae  Jong  Seo 
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5.  Peer-reviewed  Technical  Publications 


“Maximum  Likelihood  Methods  For  Time-resolved  Imaging  Through  Turbid  Media”, 

Eriksson,  B.,  Nowak,  R.,  in  Proc.  of  the  International  Conference  on  Image  Processing,  Atlanta,  GA,  October 
2006,  pp  641-644 

“Robust  Kernel  Regression  for  Restoration  and  Reconstruction  of  Images  from  Sparse,  Noisy  Data” 

Takeda,  H.,  S.  Farsiu,  P.  Milanfar,  in  Proc.  of  the  International  Conference  on  Image  Processing,  Atlanta,  GA, 
October  2006,  pp1257-1260 

“A  Decision-Theoretic  Approach  to  Transillumination  Imaging  in  Biological  Mediums”,  B.  Eriksson  And  R. 
Nowak,  Proceedings  of  the  IEEE  International  Symposium  on  Biomedical  Imaging,  April  2007. 

“Statistical  Detection  and  Imaging  of  Objects  Hidden  in  Turbid  Media  Using  Ballistic  Photons”,  S.  Farsiu,  J. 
Christofferson,  B.  Eriksson,  P.  Milanfar,  B.  Friedlander,  A.  Shakouri,  R.  Nowak,  Applied  Optics,  vol.  46,  no.  23, 
pp.  5805-5822,  Aug.2007. 

"Multi-Scale  Statistical  Detection  and  Ballistic  Imaging  Through  Turbid  Media",  S.  Farsiu  and  P.  Milanfar,  In 
proceedings  of  the  IEEE  International  Conference  on  Image  Processing  (ICIP),  San  Antonio,  TX,  Sept.  2007. 

“Resolution  Bounds  and  Reconstruction  Techniques  for  Time-Resolved  Transillumination  Imaging”,  B. 
Eriksson,  and  R.  Nowak,  Submitted  to  IEEE  Transactions  on  Image  Processing,  July  2007 


6.  Interactions/Transitions 

The  PI,  Co-PIs  and  Post-docs  have  presented  17  invited  talks  in  the  course  of  this  project.  These  have  resulted 
in  excellent  technical  interactions  and  dissemination  of  our  work.  In  many  cases,  the  talks  have  motivated 
members  of  the  audience  to  look  at  the  research  problems  we  are  posing,  and  therefore  we  believe  our  work 
has  spurred  a  healthy  amount  of  new  activity  in  application  of  statistical  signal  processing  methods  to  imaging 
problems.  In  some  cases,  these  applications  have  been  in  unexpected  areas  such  as  medical  imaging. 


6.1  Invited  Presentations 

Invited  Speakers,  Peyman  Milanfar,  Ali  Shakouri,  Robert  Nowak,  Waveforms  for  Active  Sensing  Program  Kick 
Off  Presentations,  Washington  D.C.,  Sept  2005 
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Invited  Speakers,  Peyman  Milanfar,  Ali  Shakouri,  Robert  Nowak,  Benjamin  Friendlander,  Waveforms  for  Active 
Sensing  Program,  First  Program  Review  Presentations,  Portland,  OR,  March,  2005 

Invited  Speaker,  Peyman  Milanfar,  Air  Force  Office  of  Scientific  Research  TCATS  Workshop,  Tucson,  AZ, 
2006 

Invited  Speaker,  Peyman  Milanfar,  DARPA  Waveforms  for  Active  Sensing  Workshop,  Portland,  OR,  2006 

Invited  Speaker,  Peyman  Milanfar,  European  Signal  Processing  Conference,  Florence,  Italy,  2006 

Invited  Speaker,  Peyman  Milanfar,  SIAM  Conference  on  Imaging  Science,  Minneapolis,  MN,  2006 

Invited  Speaker,  Sina  Farsiu,  Department  of  Electrical  and  Computer  Engineering,  University  of  Wisconsin, 
Madison,  Wl,  May,  2006 

Invited  Speaker,  Sina  Farsiu  ,  Waveforms  for  Active  Sensing  Program,  Second  Program  Review  Presentation, 
San-Agustin,  FL,  September,  2006 

Invited  Speaker,  Peyman  Milanfar,  International  Conference  on  Image  Processing,  Atlanta,  GA,  Oct.  2006 
Invited  Speaker,  Robert  Nowak,  International  Conference  on  Image  Processing,  Atlanta,  GA,  Oct.  2006 
Invited  Speaker,  Sina  Farsiu,  Duke  University  Eye  Center,  Durham,  NC,  Feb  2007 

Invited  Speaker,  Sina  Farsiu,  Biomedical  Engineering  Department,  Duke  University,  Durham,  NC,  Feb,  2007 
Invited  Speaker,  Sina  Farsiu,  Sony  Electronics,  San  Jose,  CA,  Feb.  2007 

Invited  Speaker,  Peyman  Milanfar  Conference  on  Applied  Inverse  Problems,  Vancouver,  BC,  Canada,  June 
2007 

Invited  Speaker,  Peyman  Milanfar  ,  AFOSR  Sensing  Program  Review,  Harvard  Univ.,  Cambridge,  MA  June 
2007 

Invited  Speaker,  Peyman  Milanfar,  SENSIP  Center,  Arizona  State  University,  Tempe,  AZ,  March  2007 
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Invited  Speaker,  Peyman  Milanfar,  Center  for  Advanced  Signal  and  Image  Sciences,  Lawrence  Livermore 
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Aside  from  the  interactions  between  the  group  members,  and  the  members  of  the  other  groups  in  this  program, 
we  have  made  contacts  with  many  prominent  researchers  in  the  academic  community.  Below  we  list  a  number 
of  university  colleagues  who  have  been  consulted: 

•  Prof.  Ines  Delfino 

Biophysics  &  NanoScience  Centre 

CNISM  -  Unita  di  Ricerca  di  Viterbo 

Dipartimento  di  Scienze  Ambientali  - 

Universita  della  Tuscia,  Largo  dell'Universita,  01100  Viterbo, 

ITALY 

•  Prof.  Robert  R.  Alfano’s  Group  especially  Dr.  M.  Alrubaiee 
Department  of  Science  And  Engineering 

City  College  and  Graduate  School,  City  University  of  New  York 
USA 

•  Prof.  Stefan  Dilhaire 

Centre  de  Physique  Moleculaire  Optique  et  Hertzienne  -  UMR  5798 
University  of  Bordeaux, 

France 
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We  exploit  recent  advances  in  active  high-resolution  imaging  through  scattering  media  with  ballistic 
photons.  We  derive  the  fundamental  limits  on  the  accuracy  of  the  estimated  parameters  of  a  mathematical 
model  that  describes  such  an  imaging  scenario  and  compare  the  performance  of  ballistic  and  conventional 
imaging  systems.  This  model  is  later  used  to  derive  optimal  single-pixel  statistical  tests  for  detecting  objects 
hidden  in  turbid  media.  To  improve  the  detection  rate  of  the  aforementioned  single-pixel  detectors,  we 
develop  a  multiscale  algorithm  based  on  the  generalized  likelihood  ratio  test  framework.  Moreover,  con¬ 
sidering  the  effect  of  diffraction,  we  derive  a  lower  bound  on  the  achievable  spatial  resolution  of  the  proposed 
imaging  systems.  Furthermore,  we  present  the  first  experimental  ballistic  scanner  that  directly  takes 
advantage  of  novel  adaptive  sampling  and  reconstruction  techniques.  ©  2007  Optical  Society  of  America 
OCIS  codes:  100.0100,  030.6600,  140.0140,  320.0320. 


1.  Introduction 

High-resolution  imaging  and  detection  of  objects  hid¬ 
den  in  a  turbid  (scattering)  medium  have  long  been 
challenging  and  important  problems  with  many 
industrial,  military,  and  medical  applications.  Al¬ 
though  turbid  media  such  as  fog,  smoke,  haze,  or 
body  tissue  are  virtually  transparent  to  radar  range 
electromagnetic  waves,  the  resolution  of  radar-based 
imaging  systems  is  often  insufficient  for  many  prac¬ 
tical  applications.  Moreover,  in  some  instances  the 
transparency  characteristics  of  certain  objects  (tar¬ 
gets)  and  the  medium  are  very  close  in  the  radar 
range  spectrum,  making  them  practically  indistin¬ 
guishable  from  each  other.  On  the  other  hand,  al¬ 
though  the  resolution  of  imaging  systems  using 
ultrashort  wavelengths  (e.g.,  x  rays)  is  desirable, 
there  exist  potential  health  hazards  for  imaging  sub¬ 
jects  and  technicians  alike. 

As  an  alternative,  imaging  systems  working  in  the 
optical-infrared  spectrum  range  (laser  scanners)  are 
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potentially  able  to  produce  high-resolution  images 
without  the  likely  health  hazards.  Unfortunately, 
even  a  very  thin  and  powerful  collimated  laser  beam 
quickly  diffuses  as  it  travels  in  turbid  media,  similar 
to  a  car’s  headlights  in  fog.  Therefore,  a  naive  ap¬ 
proach  to  optical  imaging  of  objects  hidden  inside  a 
turbid  medium  results  in  blurry  images  where  tar¬ 
gets  are  often  indistinguishable  from  each  other  or 
the  background. 

Fortunately,  the  advent  of  the  new  tunable  solid- 
state  lasers  and  ultrafast  optical  detectors  has  en¬ 
abled  us  to  acquire  high-quality  images  through 
turbid  media  where  the  resolution  is  only  limited  by 
diffraction.  Although  many  efficient  imaging  systems 
for  capturing  high-resolution  images  through  turbid 
media  have  been  proposed  throughout  the  years  [1] , 
in  this  paper  we  mainly  focus  on  ultrafast  time-gated 
or  coherent  imaging  systems  [2].  We  note  that  the 
proposed  methods  and  analysis  are  valid  and  appli¬ 
cable  for  a  great  range  of  imaging  systems  including 
optical  coherence  tomography  [3]  and  x-ray  imaging 
systems. 

Ultrafast  time-gated  imaging  is  based  on  scanning 
the  region  of  interest  (ROI)  point  by  point  by  sending 
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fast  bursts  of  optical  energy  (laser  pulses)  and  detect¬ 
ing  the  unscattered  (coherent)  photons  that  have 
passed  through  the  medium  or  reflected  from  the 
object.  Although  most  of  the  photons  in  a  laser  pulse 
are  either  randomly  scattered  (losing  their  coher¬ 
ence)  or  absorbed  as  they  travel  through  turbid  me¬ 
dia,  across  short  distances,  a  few  photons  keep  their 
coherence  and  pass  through  in  straight  lines  without 
being  scattered.  These  coherent  photons  are  com¬ 
monly  referred  to  as  the  ballistic  photons.  Aside  from 
the  diffusive  and  ballistic  photons,  the  photons  that 
are  slightly  scattered,  retaining  some  degree  of  co¬ 
herence,  are  referred  to  as  snake  photons. 

In  what  follows  in  this  paper,  we  focus  on  studying 
and  improving  the  performance  of  ballistic  imaging 
systems.  In  Section  2,  we  describe  a  statistical  model 
for  the  signal  and  noise  in  a  typical  ballistic  imaging 
scenario.  Furthermore,  we  describe  optimal  methods 
for  characterizing  the  optical  properties  of  the  scat¬ 
tering  medium  and  the  semitransparent  objects  in¬ 
side  it.  In  Section  3,  we  study  the  performance  limits 
of  optimal  single-pixel  detection  systems.  Moreover, 
we  show  that  better  detection  rates  are  achievable 
using  a  multipixel  detection  technique  based  on  the 
generalized  likelihood  ratio  test  (GLRT)  principle. 
The  effect  of  diffraction  on  the  detection  rate  is  dis¬ 
cussed  in  Section  4.  In  Section  5,  we  describe  a  lab¬ 
oratory  setup  for  detecting  ballistic  photons  and 
capturing  high-resolution  images  through  turbid  me¬ 
dia,  where  real  experimental  data  are  presented  to 
further  clarify  the  concept  of  ballistic  imaging.  In 
Subsection  5.B,  we  describe  an  adaptive  sampling 
scheme  that  effectively  reduces  the  image  acquisition 
time,  making  ballistic  imaging  more  suitable  for 
practical  applications.  A  summary  and  future  work 
directions  are  given  in  Section  6,  which  concludes  this 
paper. 

2.  Statistical  Model  for  Ballistic  Imaging  Systems 

To  have  a  better  understanding  of  the  practical  issues 
involved  in  photon-limited  imaging  via  ballistic  sys¬ 
tems,  let  us  consider  the  imaging  system  described  by 
Zevallos  et  al.  [4]  where  the  pumped  Ti:sapphire  laser 
radiates  800  nm  pulses  at  a  repetition  rate  of  1  kHz 
and  an  average  power  of  60  mW.  It  is  easy  to  show 
that  the  energy  delivered  by  the  laser  during  each 
pulse  is 

60X  10  3X  Is  ^ 

^-pulse  1000  6  x  10  J, 

and  the  energy  of  each  photon  is  computed  as 


e  =  hf=  ^  =  2.4830  X  10  19  J, 


where  h  =  6.626  X  1CT34  is  Planck’s  constant,  c  = 
299,  792,  458  m/s  is  the  speed  of  light,  and  A  = 
800  nm  is  the  wavelength.  Now  the  number  of 
photons  in  each  packet  of  energy  (pulse)  is  easily 


computed  as 


Io  = 


6  x  10 


2.4830  x  10 


=  2.4164  X  1014  photons. 


(1) 


Because  of  the  statistical  nature  of  pulse  propaga¬ 
tion,  as  a  laser  beam  travels  through  a  diffusive  me¬ 
dium,  it  is  possible  that  some  of  the  photons  emerge 
without  being  scattered.  By  selecting  these  unscat¬ 
tered  ballistic  photons  and  rejecting  the  scattered 
(diffused)  ones,  it  is  possible  to  obtain  nonblurred 
images  that  are  the  sharp  shadows  of  targets  buried 
in  the  diffusive  medium. 

Since  the  diffusive  and  ballistic  photons  have  dif¬ 
ferent  path  lengths,  a  femtosecond  laser  pulse  gen¬ 
erator  and  an  ultrafast  time  gate  can  be  paired  to 
separate  the  relatively  slow  (delayed)  diffusive  pho¬ 
tons  from  the  ballistic  ones.  We  will  say  more  on  a 
practical  setup  of  a  ballistic  photon  imaging  system  in 
Section  5.  In  what  follows  in  this  section,  we  focus  on 
modeling  the  detected  ballistic  photons  and  noise 
from  a  statistical  point  of  view. 

A.  Modeling  Received  Signal  Power 

As  expected,  in  relatively  long  distances,  the  number 
of  detected  ballistic  photons  is  extremely  small.  In¬ 
deed,  Beer’s  law  [5]  dictates  an  exponential  relation¬ 
ship  between  the  intensity  of  the  transmitted  light 
and  that  of  the  ballistic  component  as 


h  =  h  exp  |  -  L 


cL 


(2) 


In  this  expression,  I0  is  the  number  of  the  generated 
photons  in  one  laser  pulse  before  entering  the  turbid 
medium,  Ib  is  the  number  of  the  ballistic  photons  that 
survive  traveling  through  the  medium,  d  is  the  dis¬ 
tance  traveled  through  the  medium,  L  =  1/p,,  is  the 
mean  free  path  (MFP)  length  (average  distance  pho¬ 
tons  travel  before  being  scattered),  and  p,  =  ps 
+  pa  is  the  medium  extinction  factor  (the  summation 
of  scattering  and  absorptive  coefficients,  respec¬ 
tively).  From  Eqs.  (1)  and  (2),  it  is  clear  that  for  the 
laboratory  imaging  systems  with  laser  power  of  the 
order  of  the  one  described  by  Zevallos  et  al.  [4],  it  is 
fairly  unlikely  that  any  ballistic  photon  survives  im¬ 
aging  scenarios  where  the  ratio  of  d/L  is  larger  than 
—30  MFPs.  In  Appendix  A,  we  have  included  a  de¬ 
tailed  decision-theoretic  study  for  defining  the  critical 
distance  after  which  the  conventional  (non-time- 
gated)  imaging  systems  are  preferred  to  the  time¬ 
gated  ballistic  systems. 

The  exponential  drop  in  the  number  of  received 
photons  is  the  main  prohibitive  factor  for  using  such 
high-resolution  optical  imaging  systems  across  long 
distances.  In  such  imaging  scenarios,  we  are  forced  to 
rely  on  the  less  immediately  informative  (due  to  the 
inherently  severe  blur)  snake  and  diffusive  photons. 
In  recent  literature  [6,7],  an  accurate  yet  computa¬ 
tionally  manageable  mathematical  model  for  diffu- 
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sive  light  propagation  in  turbid  media  is  presented. 
Cai  et  al.  [8]  analyzed  and  experimented  on  such  an 
imaging  modality  and  Das  et  al.  [9]  and  Gibson  et  al. 
[10]  presented  some  excellent  literature  surveys  on 
the  subject  of  diffusive  imaging  systems.  However, 
imaging  systems  that  are  able  to  time  resolve  both 
ballistic  and  diffusive  photons  are  rather  expensive 
(e.g.,  a  gated  optical  intensifier  camera  costs  about 
$100,000)  and  are  not  discussed  in  this  paper.  Here, 
we  focus  on  and  derive  fundamental  performance  lim¬ 
its  for  imaging  systems  that  detect  ballistic  photons 
only.  We  exploit  these  statistical  studies  to  improve  the 
performance  of  ballistic  imaging  systems  even  in  long 
distances  where  the  signal  power  is  weak. 

It  is  important  to  note  that  because  of  the  stochas¬ 
tic  nature  of  photon  propagation,  Ib,  calculated  in  Eq. 
(2),  is  merely  the  expected  value  of  a  Poisson  random 
variable  that  estimates  the  number  of  surviving  bal¬ 
listic  photons.  Moreover,  we  assume  that  the  received 
signal  at  the  detector  is  contaminated  with  some 
amount  of  independent  Poisson  noise  due  to  shot 
noise  and  other  degrading  effects.  Therefore,  since 
the  received  signal  at  the  detector  is  the  unweighted 
summation  of  two  Poisson  random  variables,  it  can  be 
modeled  as  a  Poisson  random  process  with  the  fol¬ 
lowing  expected  value: 

I  =  I0  exp(  - 1 x,d)  +Xe=Xs+Xe, 

where  Xe  and  X,  are  the  expected  values  of  the  noise 
and  signal,  respectively.  Note  that  weighted  summa¬ 
tion  of  Poisson  random  variables  in  general  is  not 
Poissonian,  which  in  some  cases  can  be  approximated 
as  a  truncated  Gaussian  distribution  [11].  However, 
summation  of  Poisson  random  variables  with  integer 
weights  is  yet  another  Poisson  random  variable. 

B.  Characterizing  the  Optical  Properties  of  the  Medium  in 
the  Absence  of  Targets 

Accurate  characterization  of  the  scattering  medium’s 
optical  properties  is  essential  for  designing  optimal 
detectors.  Since  light  propagation  in  ballistic  imaging 
systems  is  described  by  the  single-parameter  Beer’s 
law  model,  we  are  mostly  interested  in  measuring 
(characterizing)  the  medium  or  semitransparent  ob¬ 
ject’s  extinction  factor. 

In  the  imaging  model  of  Subsection  2.A,  the  re¬ 
ceived  signal  is  modeled  as  a  Poisson  random  vari¬ 
able  with  probability  density  function 


Since  the  average  power  of  the  laser  or  the  detector 
(and  medium)  characteristics  are  assumed  not  to  be 
changing  abruptly,  to  simplify  notations,  we  assume 
that  Xei=  Xe2  =  ...  =  X„N  =  Xe  and  XSi  =  X  =  ... 
=  XSf/  =  Xs  (extension  to  the  more  general  time- 
varying  signal  and  noise  case  is  straight  forward). 
The  maximum  likelihood  (ML)  estimate  of  the 
medium’s  extinction  factor  is  given  by 


In 


NI0 


'  1°g[/’(ypC+^)] 

dfk 


=  0  =>  pf  = 


NXe-^yk 

£=1  I 


d 


Study  of  the  Fisher  information  matrix  (FIM)  de¬ 
termines  the  accuracy  of  the  above  estimation 
scheme.  Each  element  of  this  matrix  can  be  computed 
[12]  as 


a2log[/-(y|^+^)]  w  i  3XSkdXSk\ 

iJ  a0,f)0;  h\xe+xSk  (t©,:  d&J’ 

where  E  is  the  expected  value  operator  and  ®k  is  the 
Mh  parameter  of  the  model.  For  the  case  of  charac¬ 
terizing  the  extinction  factor  of  the  medium,  the  FIM 
has  only  one  element: 


NI02d2e~2lL‘d 
~  Xe+NIde-'L‘d' 

Note  that  an  unbiased  estimator  can  be  found  that 
attains  the  Cramer-Rao  bound  (CRB),  which  defines 
a  lower  bound  on  the  covariance  of  any  unbiased 
estimator  [13],  if  and  only  if  the  estimator  is  a  linear 
transformation  of  the  gradient  of  the  log-likelihood 
(score)  function  [13,14] 

log\ f(y\Ine  V‘"  \-  XM  ? 

- ,n - =  -  M*)- 


Now,  since 


dlog[f(y\I0e^  +Xe)] 


I0de^td{Xe+I^) 


1 0de  M  n 
N  Av* 


,n- 


{XH+X>:d(Xek  +XHy>-‘ 

yj 


(3) 


where  yk  is  the  Mh  measurement,  y  =  [>q,  y2,  .  .  ■  , 
Jk,  ■  ■  ■  ,  VnT,  =  K,,  Xe2,  ...  ,X„k,  ...  ,  XeNf,  and 
a;  =  [X,,  XS2,  .  .  . ,  XSk,  .  .  . ,  x,vf.  Note  that  the  laser 
emits  thousands  of  pulses  per  second  and  in  practical 
implementation  each  spatial  position  is  measured  N 
times  to  improve  the  quality  of  estimation,  and  there¬ 
fore  the  model  in  Eq.  (3)  is  presented  in  vector  form. 


it  is  clear  that  no  efficient  estimate  of  the  extinction 
parameter  can  be  found  and  such  estimates  will  al¬ 
ways  be  biased.  This  suggests  that,  in  general,  the 
lower  bound  on  the  variance  of  such  an  estimator 
cannot  be  computed  by  simply  inverting  the  Fisher 
matrix  element.  Fortunately,  we  can  numerically 
show  that  for  the  turbid  media  that  are  of  most  in¬ 
terest  to  us  [such  as  [15]  heavy  fog  (pf  =  12.5  1 
m  ’),  light  fog  (p,  =  125  1  m  '),  and  haze  (p,  = 
505.05  1  m-1)],  the  bias  component  relative  to  the 
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Fig.  1.  (Color  online)  Optimal  distance  for  calibrating  the  medium  extinction  factor  for  heavy  fog,  light  fog,  and  haze,  (a)  Experimental 
setup,  where  the  detector  is  moved  to  different  locations  [marked  by  lighter  (red)  dots]  inside  the  turbid  medium,  (b)  Bias  of  estimation 
that  is  calculated  over  60,000  Monte  Carlo  simulations,  (c)  Summation  of  squared  bias  and  variance  (solid  curves)  that  is  dominated  by 
the  variance  component  and  perfectly  fits  the  predicted  results  from  CRB  formulation  (dotted  curves)  in  short  distances. 


variance  is  small  and  can  be  ignored.  Therefore  the 
CRB  on  the  variance  can  be  expressed  as 

Xe  +/0e~M',d 

(4> 

Aside  from  theoretical  analysis,  in  practice,  this  sim¬ 
ple  closed-form  expression  of  the  lower  bound  to  the 
variance  of  the  estimate  can  help  us  design  optimal 
experiments  to  characterize  the  optical  properties  of 
the  medium  and  the  target. 

For  example,  the  CRB  analysis  helps  us  find  the 
optimal  distance  between  the  laser  and  the  detector 
for  estimating  the  medium  extinction  factor.  Figure 
1(a)  shows  the  setup  of  this  numerical  experiment, 
where  the  black  dot  represents  the  position  of  the 
laser  and  the  lighter  (red)  dots  represent  the  possible 
locations  of  the  detector.  The  optimal  distance  mini¬ 
mizing  the  lower  bound  on  the  estimator  variance  can 
be  easily  calculated  by  differentiation  of  Eq.  (4)  with 
respect  to  the  distance  (d).  Figure  1(b)  shows  the 
estimated  bias  for  this  experiment  (via  60,000  Monte 
Carlo  experiments),  which  are  small  and  negligible. 
In  Fig.  1(c),  we  have  plotted  the  summation  of  the 
numerically  experimented  bias  (squared)  and  the 
minimum  variance  (solid  curves)  and  the  CRB  (dot¬ 
ted  curves)  predicted  from  Eq.  (4),  which  perfectly  fit 
the  numerically  experimented  results  in  shorter  dis¬ 
tances.  These  plots  suggest  that,  for  calibrating 
heavy  fog,  the  optimal  distance  between  the  laser  and 
the  detector  is  less  than  100  m,  whereas  such  a  dis¬ 
tance  for  light  fog  is  of  the  order  of  a  few  hundred 
meters  and  for  haze  is  of  the  order  of  1  km.  Note  that 
the  dotted  curves  (numerically  experimented  results) 
in  Figs.  1(b) — 1(c)  are  discontinued  after  certain  dis¬ 
tances.  The  reason  for  such  discontinuity  is  that  in 
long  distances,  where  the  signal  power  is  about  the 
same  as  the  noise  level,  the  estimated  bias  is  not 
negligible  and  abruptly  tends  to  infinity.  Therefore, 
the  proposed  CRB  formulation  (4),  depending  on  the 
scattering  properties  of  the  medium,  is  only  valid  up 
to  some  distance  as  plotted  in  Fig.  1.  Practically,  this 
is  of  no  concern,  since  these  distances  are  away  from 
the  optimal  calibration  distance. 


C.  Joint  Characterization  of  the  Medium  and  the  Target’s 
Optical  Properties 

A  related  and  more  practical  problem,  namely,  char¬ 
acterizing  the  optical  properties  of  an  object  located 
inside  an  unknown  turbid  medium,  requires  two  in¬ 
dependent  sets  of  experiments.  The  first  set  of  exper¬ 
iments  is  performed  in  the  absence  of  the  object  (and 
repeated  Nt  times  to  improve  the  accuracy)  and  the 
second  set  of  experiments  is  performed  in  the  pres¬ 
ence  of  the  presumed  object  (and  repeated  N2  times). 
Figure  2  illustrates  such  an  imaging  scenario,  for 
which  we  can  easily  derive  the  ML  estimates  of  the 
medium  and  the  object  (inclusion)  extinction  factors  as 
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respectively,  where  dmc  is  the  thickness  of  the  ob¬ 
ject. 

The  general  FIM  formulation  of  Eq.  (4)  can  be  ex¬ 
ploited  for  both  of  these  imaging  scenarios.  In  this 


Fig.  2.  (Color  online)  Experimental  setup  for  characterizing  the 
optical  properties  of  the  medium  (pd  and  a  semitransparent  object 
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case,  the  CRBs  are  derived  from  the  inverse  of  a 
[2  X  2]  FIM,  the  diagonal  elements  of  which  define 
the  variance  bounds: 
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As  an  illustrative  example,  we  fixed  N1  and  N2  to  50 
each,  Xe  =  20,  and  assumed  that  semitransparent 
objects  with  extinction  factors  of  a,  =  0.124,  ix, 

=  1.24,  and  |x*.  =  12.4  and  1  m  thickness  are  present 
inside  heavy  fog.  In  Fig.  3,  we  compared  the  numer¬ 
ically  experimented  squared  bias  and  variance  (via 
5000  Monte  Carlo  simulations)  to  the  CRB  limit,  as¬ 
suming  that  the  distance  between  the  laser  and  the 
detector  are  variant  between  50  and  300  m.  The  re¬ 
sults  basically  show  that  the  numerically  experi¬ 
mented  and  CRB  values  of  the  medium  extinction 
factor  in  all  cases  are  indistinguishably  close  to 
each  other.  On  the  other  hand,  as  the  inclusive 
object  becomes  more  opaque,  the  theoretic  CRB  and 
numerically  experimented  variance  diverge  from 
each  other. 


3.  Performance  Analysis  of  Pixelwise  Optimal 
Detectors 

In  this  section,  assuming  that  the  laser,  target,  and 
turbid  medium  are  accurately  calibrated,  we  study 
the  performance  bounds  of  optimal  detectors  in  the 
presence  of  opaque  or  semitransparent  objects. 

A.  Detecting  Opaque  Objects 

In  this  subsection,  we  study  the  performance  of  the 
Neyman-Pearson  (NP)  type  statistical  test  [16]  for 
detecting  opaque  objects  hidden  in  a  turbid  medium 
versus  distance.  In  this  test,  we  basically  compare  the 
likelihood  of  the  following  two  scenarios: 

•  H0:  An  opaque  object  is  hidden  in  the  scattering 
medium,  blocking  the  laser  pulse  (i.e.,  measurements 
contain  only  noise). 

•  Hp  No  opaque  object  exists  in  the  propagation 
line  of  the  laser  pulse  (i.e.,  measurements  contain 
noise  plus  an  attenuated  laser  pulse). 
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and  therefore  the  NP  test  is  derived  by  comparing  the 
log-likelihood  ratio  to  a  threshold  as 
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Noting  that  2f=i Vk  is  yet  another  Poisson  process, 
the  probabilities  of  false  alarm  (PFA)  and  detection 
( PD )  are  computed  as 
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where  CDF  is  the  cumulative  distribution  function  of 
a  Poisson  random  variable.  Note  that,  in  some  scien¬ 
tific  communities,  the  false  alarm  is  commonly  re¬ 
ferred  to  as  a  false  positive  and  detection  is  referred 
to  as  a  true  positive. 

Figure  4(a)  shows  the  receiver  operating  character¬ 
istics  (ROC)  ( PD  versus  PFA)  curves  for  detecting 
opaque  objects  in  heavy  fog,  considering  a  detector 
with  A,.  =  20  and  a  laser  power  as  in  Eq.  (1).  This 
experiment  shows  that,  by  using  only  ballistic  pho¬ 
tons,  it  is  possible  to  reliably  detect  the  existence  (or 
absence)  of  opaque  objects  in  this  scattering  me¬ 
dium  up  to  a  distance  of  —30  MFPs.  Figure  4(b) 
shows  the  system  performance  curves  by  fixing  the 
false  alarm  rate  (PFA)  at  0.0015,  0.015,  and  0.15  val¬ 
ues  and  plotting  the  detection  rate  versus  distance 
(PD  versus  d). 

B.  Detecting  Semitransparent  Objects 

Detection  of  semitransparent  objects  is  based  on  dif¬ 
ferentiating  between  the  following  two  imaging  sce¬ 
narios: 


The  probability  density  function  of  these  two  see-  •  H0:  A  semitransparent  object  is  hidden  in  the 

narios  when  such  tests  are  repeated  N  times  are  scattering  medium,  partially  blocking  the  laser  pulse 
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Fig.  3.  (Color  online)  Comparison  of  the  bias  and  variance  from  5000  Monte  Carlo  simulations  (numerically  experimented)  and  the 
estimated  CRB  values  of  the  medium  and  the  semitransparent  target’s  optical  properties  versus  distance.  The  bias,  variance,  and  CRB 
of  the  medium  extinction  factor  are  compared  in  (a),  (b),  and  (c).  The  bias,  variance,  and  CRB  of  the  target’s  extinction  factor  are  compared 
in  (d),  (e),  and  (f). 


(i.e.,  the  measurement  is  noise  plus  signal  attenuated 
by  both  the  medium  and  the  target). 

•  IHl!:  No  semitransparent  object  exists  in  the 
propagation  line  of  the  laser  pulse  in  the  scattering 
(i.e.,  measurement  is  noise  plus  signal  attenuated  by 
the  medium). 


Xs.  —  I0e 


'M M -t(d  dine) 


where  p,s.nc  and  dinc  are  the  extinction  factor  and  the 
thickness  of  the  object,  respectively.  Based  on  this 
model,  a  NP  detection  rule  is  derived  as 


The  number  of  ballistic  photons  in  the  attenuated 
signal  that  travel  through  both  the  medium  and  the 
semitransparent  object  is  calculated  as 
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Fig.  4.  (Color  online)  (a)  ROC  plots  at  different  distances  for  detecting  opaque  objects  in  heavy  fog  (p,  =  12.5  1  m  x)  and-X),  =  20.  (b)  By 
fixing  the  P FA  at  different  values,  the  detection  rate  (PD)  is  plotted  versus  the  distance. 


The  probabilities  of  false  alarm  and  detection  are 
computed  as 
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Figure  5  shows  the  ROC  curves  for  detecting  a 
semitransparent  object  in  heavy  fog.  In  this  experi¬ 
ment  using  a  laser  and  detector  similar  to  the  ones  in 
Subsection  3.A,  the  distance  was  fixed  at  300  m, 
which  according  to  Fig.  4  delivers  almost  perfect  de¬ 
tection  for  opaque  objects.  Figure  5(b)  shows  the  sys¬ 
tem  performance  curves  by  fixing  the  false  alarm  rate 
(PFa)  at  0.00015,  0.0015,  and  0.015  values  and  plot¬ 
ting  the  detection  rate  versus  the  object’s  extinction 
factor  (PD  versus  As  expected,  this  experiment 
shows  that  the  detection  performance  deteriorates  as 
the  object  becomes  less  opaque. 

C.  Multipixel  GLRT  Detection 

As  explained  in  Subsection  2. A,  in  ballistic  imaging 
the  field  of  view  is  scanned  at  multiple  points  to  cre¬ 
ate  a  2-D  image  of  the  objects  in  the  ROI.  In  this 
subsection,  we  propose  an  effective  algorithm  that 


(a)  (b) 

Fig.  5.  (Color  online)  (a)  ROC  plots  for  detecting  transilluminative  objects  at  300  m  distance  in  heavy  fog  (p,  =  12.5-1  m)  and 
Xe  =  20.  (b)  By  fixing  the  PFA  at  different  values,  detection  rate  (PD)  is  plotted  versus  the  object’s  transparency  (p,im). 
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exploits  the  spatial  correlation  of  the  nearby  samples 
in  a  multipixel  imaging  scenario  to  improve  on  the 
performance  of  the  single-pixel  optimal  detectors  de¬ 
veloped  in  the  previous  section. 

The  proposed  multipixel  detection  technique  gen¬ 
eralizes  the  single-pixel  detection  techniques  and 
performs  tests  on  superpixels,  which  are  the  collec¬ 
tive  intensities  of  a  set  of  neighboring  pixels  in  size 
and  shape  of  the  hidden  objects.  However,  since  in 
general  the  size  and  shape  of  the  hidden  objects  is  not 
known  a  priori,  we  develop  a  GLRT-based  algorithm 
that  simultaneously  tests  the  existence  and  also  es¬ 
timates  the  shape  and  size  of  the  objects  hidden  in 
turbid  media. 

The  outline  of  the  proposed  GLRT  algorithm  is 
illustrated  by  an  example  in  Fig.  6.  First,  for  a  given 
(fixed)  false  alarm  rate  the  optimal  detectors  devel¬ 
oped  in  the  previous  section  are  exploited  to  test  the 
existence  or  absence  of  objects  at  each  individual 
pixel.  As  an  illustrative  example,  this  test  is  applied 
to  the  central  pixel  (shaded)  of  Fig.  6(a),  where  the 
measured  pixel  value  (0.4)  is  compared  with  the  NP 


-0.1 

(a)  (b) 


1.9 

/ 

(c)  (d) 


Fig.  6.  (Color  online)  Illustrative  example  showing  the  outline  of 
the  proposed  multiscale  GLRT  algorithm.  The  check-marked  sec¬ 
ond  scale  gives  the  highest  confidence  value  for  the  central  pixel. 


test  threshold  (0.5).  Of  course,  the  greater  the  dis¬ 
tance  of  the  measurement  from  the  threshold,  the 
more  confident  we  are  in  the  accuracy  of  the  test 
result.  Next,  we  integrate  the  gray-level  values  of  all 
immediate  neighboring  pixels,  and  in  effect  consider 
them  as  one  superpixel,  as  illustrated  in  Fig.  6(b). 
Since  the  false  alarm  rate  is  fixed  for  all  scales,  the 
decision  threshold  is  different  than  the  threshold  cal¬ 
culated  in  the  previous  step,  which  is  recalculated 
based  on  the  gray-level  value  of  the  superpixel.  In  the 
next  steps,  we  repeat  this  process  by  fixing  the  false 
alarm  rate  and  considering  larger  neighborhoods. 
The  generalized  NP  test  for  these  steps  is  formulated 
as  follows: 


scale 

J m,l 


\og(rale)+Nscaljcs 


(13) 


where  ym/cale  is  the  summation  of  the  pixel  values  in 
the  Nscate  =  N( 2  X  scale  -  l)2  pixels  neighborhood 
around  the  pixel  at  position  [m,  l\.  Our  confidence  in 
the  decision  made  on  each  scale  is  simply  defined  as 
the  distance  between  the  summation  of  measure¬ 
ments  in  the  superpixel  and  that  of  the  threshold  set 
by  the  GLRT: 
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Finally,  we  decide  on  the  presence  or  absence  of  the 
object  at  a  particular  pixel  based  on  the  test  result  of 
the  scale  that  shows  the  highest  confidence  value. 
Note  that  the  optimal  scale  is  not  unique  for  all  pix¬ 
els,  as  finer  scales  are  more  suitable  for  pixels  located 
on  the  texture  or  edge  areas,  and  coarser  scales  are 
more  suitable  for  the  pixels  located  in  flat  areas.  The 
memory  requirements  of  this  technique  are  indepen¬ 
dent  of  the  maximum  scale  number,  since  we  only 
need  to  keep  the  original  image,  the  last  estimated 
image,  and  the  corresponding  confidence  values. 

To  have  a  better  understanding  of  the  proposed 
multiscale  GLRT  technique  and  its  performance,  we 
set  up  an  illustrative  controlled  imaging  scenario. 
Figure  7(a)  shows  an  ideal  (noiseless  and  determin¬ 
istic)  image  of  objects  of  different  sizes  and  shapes.  To 
depict  an  experiment  at  the  limit  distance  where  the 
signal  of  interest  is  weak,  we  consider  an  imaging 
scenario  in  which  the  average  number  of  received 
ballistic  photons  for  each  pixel  is  one  photon.  Figure 
7(b)  shows  such  Poisson  random  signals  (free  of  ad¬ 
ditive  noise  effect). 

Detection  of  such  signals  becomes  more  difficult 
when  we  consider  the  system  noise  as  illustrated  in 
Figs.  8(a)  and  9(a),  where  the  Poisson  noise  variances 
(mean)  are  20  and  40,  respectively.  Figs.  8(b)  and  9(b) 
are  images  reconstructed  by  implementing  the  point- 
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(a)  (b) 

Fig.  7.  (a)  Ideal  deterministic  and  noise-free  image  of  four  objects 
of  different  sizes  and  shapes,  (b)  Corresponding  image  as  a  Pois- 
sonian  noise-free  stochastic  signal,  with  Xs  =  1. 


by-point  single-pixel  detection  techniques,  consider¬ 
ing  a  false  alarm  rate  of  0.00125,  where  none  of  the 
objects  are  correctly  identified.  On  the  other  hand, 
Figs.  8(c)  and  9(c)  are  the  results  of  exploiting  the 
multiscale  GLRT  techniques,  showing  a  considerably 
more  accurate  detection  of  such  objects.  Figures  8(d) 
and  9(d)  illustrate  the  scale  from  which  each  pixel  in 
the  final  images  of  Figs.  8(c)  and  9(c),  respectively, 
are  selected.  Note  that,  as  expected,  the  pixels  in  the 
flat  area  are  selected  from  the  coarser  scales,  whereas 
the  pixels  on  the  edge  areas  are  selected  from  the 
finer  scales.  Figures  8(e)  and  9(e),  show  the  confi¬ 
dence  in  the  detection  result  (14)  with  respect  to  the 
corresponding  pixels.  These  figures  show  higher  con¬ 
fidence  levels  in  the  flat  and  less  confidence  in  the 


edge  areas.  Also,  in  Fig.  9(e)  we  see  that  the  area  with 
the  lowest  confidence  is  the  place  where  most  mis- 
classifications  happen.  This  is  good  news,  since  to 
increase  the  detection  rate,  we  may  opt  to  do  a  second 
(and  faster)  round  of  scans,  sampling  only  on  these 
very  low-confidence  regions.  In  Figs.  8(f)  and  9(f),  we 
plot  the  misclassification  rates  at  each  scale  (curve), 
and  compare  it  with  the  overall  multiscale  rate  (line). 
These  numerically  experimented  plots  show  that  the 
performance  of  the  proposed  pixelwise  GLRT  tech¬ 
nique  (depending  on  the  noise  level)  is  either  very 
close  to  [Fig.  8(f)]  or  even  better  than  [Fig.  9(f)]  the 
best  fixed-scale  technique.  In  Figs.  10(a)  and  10(b), 
the  performance  of  the  single-pixel  detection  tech¬ 
nique  is  compared  with  the  multiscale  ones  via  their 
corresponding  ROC  curves.  Once  again,  the  multi¬ 
scale  technique  shows  the  best  or  close  to  the  best 
performance. 

4.  Diffraction  Effects 

So  far  in  this  paper,  all  detection  tests  and  related 
performance  analysis  were  derived  based  on  a  sim¬ 
plified  model  of  light  propagation  that  ignores  diffrac¬ 
tion.  Although  such  approximation  works  well  for 
many  practical  applications,  it  is  not  a  suitable  model 
for  detecting  or  imaging  relatively  small  sized  objects. 
In  this  section,  we  present  statistical  analysis  of  the 
resolution  limits  in  ballistic  imaging  systems  by  de¬ 
fining  the  smallest  size  of  resolvable  objects  in  a  tur¬ 
bid  medium  at  given  false  alarm  and  detection  rates. 


2  4  6  8  10  12  14  16  18  20  22 
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Fig.  8.  (Color  online)  Application  of  the  proposed  multiscale  GLRT  technique  for  improving  the  detection  rate,  (a)  Result  of  adding  Poisson 
noise  ( Xe  =  20)  to  Fig.  7(a).  (b)  Result  of  the  single-pixel  detection,  (c)  Result  of  the  proposed  multiscale  detection  technique,  (d)  Image  that 
corresponds  to  the  selected  scales  for  the  image  shown  in  (c).  (e)  Corresponding  confidence  values,  (f)  Misclassification  probability  in 
different  scales. 
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Fig.  9.  (Color  online)  Application  of  the  proposed  multiscale  GLRT  technique  for  improving  the  detection  rate,  (a)  Result  of  adding  Poisson 
noise  (Xe  =  40)  to  Fig.  7(  a  ).  (b)  Result  of  the  single-pixel  detection,  (c)  Result  of  the  proposed  multiscale  detection  technique,  (d)  Image  that 
corresponds  to  the  selected  scales  for  the  image  shown  in  (c).  (e)  Corresponding  confidence  values,  (f)  Misclassification  probability  in 
different  scales. 


Our  study  is  the  continuation  of  previous  work  [17], 
generalized  by  considering  the  effects  of  the  turbid 
medium  and  ballistic  imaging  setup. 

Figure  11  is  a  1-D  illustration  of  the  diffraction, 
which  is  described  as  “any  deviation  of  light  rays  from 
rectilinear  paths  which  can  not  be  interpreted  as  re¬ 
flection  or  refraction”  [18].  In  this  figure,  the  dashed 
(red)  line  represents  the  case  in  which  light  propa¬ 
gates  in  a  straight  line  creating  a  sharp-edged 


shadow  of  a  hidden  opaque  object  at  the  detector.  The 
solid  (green)  curve,  on  the  other  hand,  illustrates  a 
more  realistic  case  in  which  the  object’s  shadow  ap¬ 
pears  blurry  at  the  detector  as  a  result  of  diffraction. 

The  blur  induced  by  diffraction  can  be  calculated 
from  the  Helmholtz-Kirchhoff  wave  propagation 
equations  [19].  However,  for  the  experimental  setups 
that  are  of  most  interest  to  us  with  respect  to  the 
scope  of  this  paper,  the  diffraction  effect  can  be  taken 


Fig.  10.  (Color  online)  Application  of  the  proposed  multiscale  GLRT  technique  for  enhancing  the  detection  rate.  ROC  plots  for  the 
proposed  multiscale  detection  scenario  in  the  imaging  scenarios  of  Figs.  8  and  9  (with  25  Monte  Carlo  experiments)  are  shown  in  (a)  and 
(b),  respectively.  The  numerical  labels  “1,  4,  .  .  .  ,  23”  correspond  to  the  scale  at  which  detection  tests  are  performed,  and  the  plot  labeled 
“Final”  represents  the  performance  of  the  proposed  multiscale  (fused)  technique. 
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Fig.  11.  (Color  online)  Shadow  of  an  opaque  object  illuminated  by 
a  homogeneous  widespread  light  beam.  The  dashed  (red)  curve 
represents  the  intensity  of  the  measured  light  ignoring  diffraction. 
The  solid  (green)  curve  represents  the  diffraction-induced  PSF. 


into  account  by  convolving  the  expected  signal  value 
with  an  appropriate  point-spread  function  (PSF)  that 
describes  the  blurring  effect  of  an  object  estimated 
from  the  Fraunhofer  approximation.  For  a  circular 
opaque  object,  such  a  PSF  is  given  by 
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where  k  =  2tt/A  is  the  wavenumber,  J  x( )  is  the  order 
1  Bessel  function  of  the  first  kind,  p  is  the  radius  of 
the  opaque  object,  r  is  the  radius  coordinate  in  the 
detector  plane,  and  z  is  the  distance  between  the 
object  and  the  detector  [20] .  Note  that  the  Fraunhofer 
approximation  is  only  valid  when  z  >>  4p2/A,  and 
therefore  in  this  paper  we  only  consider  far-field  im¬ 
aging  scenarios.  Ignoring  the  effect  of  a  turbid  me¬ 
dium  and  considering  homogeneous  illumination 
with  intensity  T  at  the  object  plane,  the  radially  sym¬ 


metric  intensity  of  the  detected  signal  at  the  radius 
coordinate  r  is  simply  given  by  I  =  I'H{r).  Figure 
12(a)  shows  the  diffraction  pattern  of  a  circular 
opaque  object  with  a  distance  ofz  =  100  m  from  the 
detector,  illuminated  by  unit-intensity  (/'  =  1)  light 
with  800  nm  wavelength. 

In  the  imaging  scenarios  considered  in  this  paper, 
the  detected  signal  is  further  attenuated  by  the  tur¬ 
bid  medium,  and  the  expected  value  of  the  signal 
intensity  at  the  radius  coordinate  r  and  distance  z 
from  the  object  plane  can  be  approximated  as 

1  =  1’  exp(-p,(z)i?(r), 


where  we  have  ignored  the  fact  that  due  to  diffraction 
some  parts  of  the  wavefront  travel  slightly  longer 
distances.  Note  that  in  practice,  due  to  the  far-field 
imaging  assumption,  such  variance  in  attenuation  is 
small.  This  effect  is  shown  in  Fig.  12(b),  where  the 
path  lengths  L1  and  L2  are  practically  equal  if  the 
distance  between  the  opaque  object  and  the  detector 
(z)  is  significantly  larger  than  the  PSF  spread.  The 
detection  problem  associated  with  the  signal  model 
defined  above  is  described  in  the  following  two  imag¬ 
ing  scenarios: 

•  H0:  An  opaque  object  of  unknown  but  small  size 
(p  >  0)  is  hidden  in  the  scattering  medium,  blocking 
and  blurring  the  laser  pulse  (i.e.,  measurements  con¬ 
tain  noise  plus  attenuated  and  blurred  laser  pulse). 

•  Hp  No  opaque  object  exists  (p  =  0)  in  the  prop¬ 
agation  line  of  the  laser  pulse  (i.e.,  measurements 
contain  noise  plus  attenuated  laser  pulse). 

The  above  GLRT  detector  is  different  than  the  NP 
detectors  of  Section  3  since  the  size  of  the  object  is 
now  assumed  to  be  unknown.  Following  Eq.  (2),  the 
expected  value  of  the  intensity  in  the  absence  of  the 
object  (Hi)  is  given  by 


I(k,  0)=Iae-'L‘d+Xe. 


Distance  (m) 


(a)  (b) 

Fig.  12.  (Color  online)  (a)  1-D  slice  of  the  diffraction  pattern  of  a  circular  object  of  radius  3  mm  at  100  m  distance  and  800  nm  wavelength, 
(b)  1-D  slice  of  the  imaging  scenario,  where  z  is  the  distance  between  the  opaque  circular  object  (radius  p)  and  the  detector,  d  is  the  distance 
between  the  laser  and  the  detector.  The  pass  length  L1  —  L2  when  z  is  very  long. 
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The  intensity  of  the  signal  in  the  presence  of  the 
object  (D-0o)  is  estimated  as 


I(k,  p)=I0e^ldH(rk)+Xe,  (16) 

where  p  is  the  estimate  of  the  opaque  object’s  radius 
and  rk  is  the  radial  distance  of  the  Mh  pixel  from  the 
axes  passing  through  the  center  of  the  object.  The 
unknown  radius  of  the  opaque  object  is  estimated  as 


w  e  l<h’">\I(k,  p) J  ' 
p  =  arg  max  ]]  — 


3V 


(17) 


The  above  ML  estimate  of  the  radius  is  solved  by 
numerical  optimization,  where  we  discretize  p  over 
an  assumed  range  of  values  p[g],  g  =  1,  G,  and 
compute  the  cost  function, 


e[S]  =  2  37  !og{7(&,  p[g\)  1  -  I(k,  p[g\).  (18) 

k  =  \ 

The  value  of  g  for  which  g\g\  takes  on  the  largest 
value  is  gmax,  and  finally  the  GLRT  detection  statis¬ 
tics  is  given  by 
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As  an  illustrative  example,  by  fixing  the  false 
alarm  rate  at  PFA  =  0.1,  the  noise  level  atX„  =  20,  and 
assuming  a  large  detector  that  detects  all  the  light, 
regardless  of  the  distance  or  size  of  the  object,  we 
used  the  above  GLRT  framework  to  search  for  the 
smallest  detectable  object  size  at  different  distances 
and  detection  rates  in  heavy  fog  (jq  =  12.5  1  m-1). 
Figure  13  illustrates  the  result  of  this  experiment, 
where  as  expected  the  size  of  detectable  objects  first 
rises  as  the  distance  increases. 


200 


Distance  (m) 

Fig.  13.  (Color  online)  Detection  rate  versus  the  (unknown) 
opaque  circular  target’s  radius  and  the  distance  between  the  laser 
and  the  detector  considering  the  diffraction  limit  with  PFA  =  0.1 
and  Xe  =  20  in  heavy  fog  (p*  =  12.5-1  m_1). 
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Fig.  14.  (Color  online)  Laser  setup  at  the  Ballistic  Imaging  Lab¬ 
oratory  at  the  University  of  California,  Santa  Cruz. 


5.  Laboratory  Setup  and  Experiments 

A.  Conventional  Ballistic  Imaging  Experimentation 

In  our  experiments,  to  generate  ultrashort  optical 
pulses  we  used  a  Coherent  Mira  900  Ti:sapphire  tun¬ 
able  femtosecond  laser  pumped  by  an  8  W  pump 
(Verdi  V-8).  At  the  output  this  laser  generates  an 
average  power  of  ~1  W  with  pulses  of 200  fs  duration, 
13  ns  repetition  period,  and  830  nm  wavelength. 

As  shown  in  Fig.  14,  each  laser  pulse  passes 
through  a  X./2  plate  and  is  incident  on  a  polarizing 
beam  splitter  that  divides  the  pulse  into  two  copies, 
one  used  for  triggering  the  ultrafast  time  gate  while 
the  other  passes  through  the  scattering  medium 
(which  is  modeled  by  two  sets  of  solid  diffusers  lo¬ 
cated  in  front  and  back  of  the  target).  Rotation  of  the 
X/2  plate  determines  the  power  ratio  between  the  two 
pulses,  and  we  experimentally  determined  that  the 
best  results  are  achieved  in  a  near  50%/50%  splitting 
ratio.  After  passing  through  the  diffusers  and  target, 
the  ballistic  photons  are  incident  on  the  gate  at  ex¬ 
actly  the  same  time  as  the  triggering  pulse  and  pass 
through  the  ultrafast  time  gate,  where  due  to  the 
phase  and  polarization  difference  the  scattered  pho¬ 
tons  are  rejected.  In  practice,  the  triggering  pulse 
timing  is  controlled  by  a  delay  line,  which  increases 
or  decreases  the  optical  path  length,  using  a 
computer-controlled  translation  stage. 

The  ultrafast  time  gate  used  is  a  nonlinear  crystal, 
(3-barium  borate  (BBO)  [21],  which  utilizes  a  two- 
photon  process  such  that  the  gating  time  can  be  as 
short  as  the  laser  pulse  width.  Additionally,  by 
slightly  changing  the  incident  angles  of  the  two 
pulses  on  the  nonlinear  crystal,  the  time-gated  result 
can  be  spatially  separated  from  the  background  sig¬ 
nal,  greatly  increasing  the  signal-to-noise  ratio.  This 
effect  is  sometimes  referred  to  as  background-free 
cross  correlation  [22] .  The  energy  of  the  ballistic  pho¬ 
tons  are  then  measured  by  a  silicon  detector  and  a 
lock-in  amplifier.  The  entire  setup  implemented  at 
the  ultrafast  imaging  laboratory  at  the  University  of 
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Fig.  15.  (Color  online)  Comparison  of  diffusive  and  ballistic  im¬ 
aging.  (a),  (b)  Two  diffusive  (no  time  gating)  scans,  (c),  (d)  Two 
corresponding  ballistic  (time-gated)  scans  through  five  solid 
ground  glass  diffusers. 


California,  Santa  Cruz,  is  controlled  using  LabView 
and  a  general  purpose  interface  bus  (GPIB)  bus. 

Figure  15  shows  the  results  of  two  imaging  exper¬ 
iments  where  the  objective  is  to  read  the  text  written 
on  transparency  sheets  by  a  ballpoint  pen,  sur¬ 
rounded  by  a  total  of  five  solid  glass  (Thorlabs  ground 
glass,  DG10-220)  diffusers.  The  thickness  of  these 
diffusers  is  2  mm  each,  with  the  MFP  of  0.73  mm  (7.3 
MFP  total).  We  also  note  that  the  dynamic  range  of 
our  system  is  approximately  100  dB.  Figures  15(a) 
and  15(b)  show  the  result  of  scans  in  the  absence  and 
Figs.  15(c)  and  15(d)  show  the  result  of  scans  in  the 
presence  of  the  ballistic  time  gate,  without  any  post¬ 
processing.  To  acquire  the  nongated  images,  the  gate 
was  removed  and  the  detector  repositioned.  Note  that 
Figs.  15(a)  and  15(c)  illustrate  raw  data  (without  any 
postprocessing)  from  the  imaging  system,  whereas 
the  results  shown  in  Figs.  15(b)  and  15(d)  are  each 
upsampled  by  a  factor  of  3  via  the  bicubic  interpola¬ 
tion  technique.  These  results  show  that,  although 
ballistic  images  are  noisier  than  the  non-time-gated 
(diffusive)  ones,  they  are  preferable  since  they  are 
virtually  blur  free. 

B.  Adaptive  Sampling  Experimentation 
As  explained  in  the  previous  sections,  in  ballistic  im¬ 
aging,  2-D  images  of  the  objects  in  the  scattering 
media  are  created  by  a  relatively  time-consuming 
point-by-point  scanning  scheme  in  which  the  field  of 
view  (FOV)  is  sampled  at  regularly  spaced  locations. 
For  instance,  in  a  typical  laboratory  setup  with  a 
mechanical  translation  stage  (Fig.  14),  creating  a 
256  X  256  image  (i.e.,  sampling  at  65,536  points) 


takes  about  4  h,  which  might  be  prohibitively  long  for 
many  real-world  applications.  Although  such  exces¬ 
sive  time  can  be  reduced  if  the  mechanical  transla¬ 
tion  stage  is  replaced  by  a  more  expensive  optical  one, 
faster  scans  are  always  desired,  and  moreover,  for 
many  applications,  the  total  number  of  pulses  deliv¬ 
ered  in  a  given  time  period  is  limited  by  the  average 
delivered  energy  due  to  health  concerns. 

By  making  some  simplifying  assumptions  about 
the  objects  of  interest  (e.g.,  piecewise  constancy),  ir¬ 
regular  scan  strategies,  such  as  sampling  sparsely  in 
the  low-frequency  areas  and  densely  in  the  high- 
frequency  (edge  or  textured)  areas,  are  shown  to  be 
useful  in  reducing  the  imaging  time.  Recently,  two 
related  techniques,  namely,  compressive  sensing 
[23,24]  and  active  learning  [25]  (adaptive  sampling) 
were  proposed  to  reduce  the  number  of  samples  re¬ 
quired  to  achieve  certain  reconstruction  accuracy 
with  respect  to  the  regular  (passive)  scanning  tech¬ 
nique.  We  note  that  it  is  the  sparsity  of  the  signal  of 
interest  (in  a  given  overcomplete  dictionary  of  bases) 
that  enables  such  techniques  to  gather  sufficient  in¬ 
formation  to  achieve  optimal  (if  not  perfect)  recon¬ 
struction  in  the  presence  of  noise,  even  when  the 
sampling  rate  is  lower  than  the  Nyquist  rate  [26].  In 
the  compressive  sensing  technique,  random  projec¬ 
tions  of  the  signal  of  interest  onto  an  overcomplete  set 
of  basis  functions  are  sequentially  recorded.  Such 
random  projections  in  practical  optical  imaging  sce¬ 
narios  can  be  implemented  by  passing  a  wide-held 
beam  through  binary  masks  with  a  random  pattern 
[in  practice,  a  digital  micromirror  device  can  be  used 
to  generate  the  random  basis  patterns  [27]].  Unfor¬ 
tunately,  in  the  ballistic  imaging  setup,  creating  a 
wide-held  beam  is  not  easy.  Moreover,  diffraction  lim¬ 
its  the  resolution  of  the  binary  mask  and  therefore 
implementing  a  compressive-sensing-based  ballistic 
imaging  system  is  not  trivial.  On  the  other  hand,  in 
the  following  we  show  that  adaptive  sampling  tech¬ 
niques  can  be  readily  exploited  for  ballistic  imaging 
purposes. 

We  have  implemented  adaptive  sampling  as  a  two- 
step  process  [28].  In  the  hrst  step,  we  regularly  sam¬ 
ple  the  FOV  space  at  N/2  points,  where  N  is  the  total 
number  of  samples  that  we  plan  to  collect.  We  use 
these  N/2  measurements  to  create  a  pilot  estimate  of 
the  unknown  FOV.  In  the  next  step,  the  remainder  of 
the  N/2  points  are  used  to  sample  the  FOV  on  the 
edge  areas  of  the  estimated  image.  It  can  be  shown 
that  the  decay  rate  of  the  mean  square  error  for  piece- 
wise  constant  images  is  0(N~1/2)  and  0(N  ')  for  the 
passive  and  active  sampling  techniques,  respectively 
[28], 

We  also  note  that  active  learning  relies  on  accurate 
adaptive  image  reconstruction  algorithms  to  recon¬ 
struct  the  unknown  images  from  the  irregular  sam¬ 
ples  of  the  FOV.  In  our  implementation,  we  used  an 
image  reconstruction  method  based  on  maximum 
a  posteriori  (MAP)  with  bilateral  total  variation  prior 
(regularizer)  [29].  The  general  formulation  of  this 
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technique  is  presented  as  follows: 


X(t)  =  arg  min 
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where  X  of  size  [ML  X  1]  is  a  vector  representing  the 
reconstructed  image  of  size  [M  X  L ]  after  lexico¬ 
graphic  ordering,  and  Z  of  size  [ML  X  1]  is  a  vector 
that  stores  the  N  <  ML  measurements.  In  this  vector, 
the  elements  that  correspond  to  those  pixels  in  X  for 
which  no  measurement  is  available  are  filled  with 
zeros.  The  matrices  Sj  and  Sym  are  the  operators 
corresponding  to  shifting  the  image  represented  byX 
by  l  pixels  in  the  horizontal  direction  and  m  pixels  in 
the  vertical  direction,  respectively.  The  scalar  weight, 
0  <  a  <  1,  is  applied  to  give  a  spatially  decaying  effect 
to  the  summation  of  the  regularization  terms,  which 
in  effect  represent  derivatives  across  multiple  reso¬ 
lution  scales.  Matrix  A  of  size  [ML  X  ML]  is  a  diag¬ 
onal  matrix  whose  values  are  chosen  in  relation  to 
our  confidence  in  the  measurements  that  contributed 
to  make  each  element  of  Z  (diagonal  elements  corre¬ 
sponding  to  pixels  for  which  no  measurement  is  avail¬ 
able  are  replaced  with  zeros).  The  regularization 
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Fig.  16.  Comparison  of  passive  and  active  imaging,  (a)  Result  of 
scanning  the  turbid  medium  on  a  regular  128  X  128  grid  (16,384 
dense  passive  sampling),  (b)  Result  of  scanning  the  turbid  medium 
on  a  regular  32  X  32  grid  (1024  sparse  passive  sampling)  followed 
by  interpolation  via  bicubic  interpolation  to  reconstruct  the  image 
on  a  256  X  256  grid,  (c)  Result  of  scanning  the  turbid  medium  on 
an  irregular  grid  (984  sparse  adaptive  sampling)  followed  by  in¬ 
terpolation  via  adaptive  interpolation  to  reconstruct  the  image  on 
a  256  X  256  grid,  (d)  Distribution  of  the  984  irregular  samples. 


parameter,  i|»,  is  a  scalar  for  properly  weighting  the 
first  term  (data  fidelity  cost)  against  the  second  term 
(regularization  cost). 

To  validate  the  applicability  of  this  adaptive  sam¬ 
pling  and  reconstruction  technique  versus  the  com¬ 
mon  passive  sampling  technique,  we  performed  the 
following  experiment.  A  metal  washer  was  imaged 
through  a  ground  glass  diffuser  (Thorlabs  ground 
glass,  DG10-220)  via  the  ballistic  imaging  setup  of 
Fig.  14.  Figure  16(a)  shows  the  result  of  scanning  the 
medium  on  a  128  X  128  (16,384  total)  regularly  sam¬ 
pled  grid.  Figure  16(b)  shows  the  same  image  sam¬ 
pled  on  a  regular  32  X  32  (1024)  grid  and  then 
upsampled  by  the  bicubic  interpolation  method  to  a 
dense  256  X  256  grid.  The  alternative  sampling 
strategy  was  performed  by  exploiting  the  same  ex¬ 
perimental  setup  (distance,  turbid  medium,  target), 
where  a  total  of  950  irregularly  sampled  data  points 
were  collected  in  the  said  two-step  adaptive  process. 
Figure  16(c)  shows  the  result  of  such  an  adaptive 
sampling  scheme  after  upsampling  to  the  256  X 
256  grid  by  the  proposed  adaptive  MAP-based  in¬ 
terpolation  method.  The  spatial  position  of  the  984 
adaptive  samples  are  marked  as  white  dots  on  a 
256  X  256  grid  in  Fig.  16(d),  which,  as  expected,  is 
considerably  denser  on  the  edge  areas. 

6.  Conclusion  and  Future  Work 

In  this  paper,  we  have  studied  a  technique  for  cap¬ 
turing  high-resolution  images  through  turbid  media. 
This  approach  was  based  on  separating  the  unscat¬ 
tered  (ballistic)  photons  from  the  diffused  ones  by 
implementing  an  ultrafast  time-gating  system.  The 
novelty  of  this  paper  is  in  combining  the  recent  ad¬ 
vances  in  optical  science  with  the  novel  image  pro¬ 
cessing  and  statistical  signal  processing  techniques. 
We  studied  the  resolution  limits  of  such  a  system  that 
were  close  to  diffraction  (Rayleigh)  limits  for  longer 
distances.  We  derived  the  fundamental  limits  on  the 
accuracy  of  the  estimated  extinction  parameters  of  an 
unknown  turbid  medium  and  the  targets  inside  it. 
This  study  also  guided  us  toward  the  most  efficient 
experiments  (with  respect  to  both  time  and  accuracy) 
for  calibrating  the  model  parameters  of  the  unknown 
turbid  medium  as  well  as  the  optical  properties  of  the 
target  (which  can  be  used  to  identify  and  categorize 
it).  Our  results  showed  that  for  a  medium  of  practical 
interest,  namely,  heavy  fog,  optical  parameters  can 
be  estimated  with  high  accuracy.  We  used  the  said 
model  to  derive  optimal  statistical  tests  for  detecting 
objects  hidden  in  turbid  media.  Performance  analysis 
was  carried  out  by  computing  ROC  curves  for  the 
proposed  optimal  tests,  showing  that,  by  considering 
only  the  ballistic  photons,  we  are  able  to  detect 
opaque  objects  hidden  in  heavy  fog  in  the  range  of 
approximately  380  m  (i.e.,  30  MFPs).  The  detection 
rate  of  the  semitransparent  objects  is  shown  to  be 
slightly  less  than  this  distance.  Also,  real  experi¬ 
ments  attested  to  the  fact  that  ballistic  imaging,  es¬ 
pecially  in  longer  distances,  is  difficult,  and  therefore 
we  developed  a  multiscale  GLRT  algorithm  to  im¬ 
prove  the  detection  rate  in  such  scenarios.  To  reduce 
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the  data  acquisition  time  that  is  essential  for  many 
real-world  applications,  we  implemented  an  adaptive 
sampling  scheme  that  significantly  reduced  the  data 
acquisition  time. 

As  for  future  work,  one  may  exploit  temporal,  spa¬ 
tial,  and  wavelength  diversity  and  coding  for  ballistic 
imaging.  We  can  study  the  array  imaging  framework 
where  multiple  emitters  will  transmit  coordinated 
pulses  of  light,  and  for  their  part,  a  collection  of 
photon-detecting  elements  will  gather  the  received 
data  and  compute  a  resultant  image.  Furthermore, 
there  is  the  possibility  of  analyzing  various  ways  of 
coding  these  pulses  (e.g.,  transmitting  sequences  of 
pulses  as  is  done  in  ultrawideband  communications). 

Moreover,  we  believe  that  detection  techniques 
that  exploit  all  ballistic,  snake,  and  diffused  photons 
[30]  (what  we  term  holistic  imaging  and  detection) 
enable  detection  of  larger  objects  at  significantly 
longer  range.  More  theoretical  and  experimental 
work  needs  to  be  done  to  design  a  (near)  optimal  yet 
practical  solution  to  this  important  problem. 

Appendix  A:  Decision-Theoretic  Resolution  Bounds 

As  explained  throughout  this  paper,  the  diffraction- 
limited  resolution  of  ballistic  imaging  systems  [e.g., 
Figs.  15(c)  and  (d)]  makes  them  appealing  for  imag¬ 
ing  in  relatively  short  distances.  However,  in  rela¬ 
tively  long  distances  the  ballistic  signal  is  too  weak 
and  we  are  bound  to  rely  on  the  blurry  but  higher 
signal-to-noise  ratio  (SNR)  images  of  conventional 
imaging  systems  [e.g.,  Figs.  15(a)  and  (b)].  In  this 
Appendix,  we  adapt  a  decision-theoretic  approach  to 
the  resolution  bounds  and  search  for  the  critical  dis¬ 
tance  after  which  the  ballistic  imaging  systems  are  of 
no  practical  advantage  compared  with  the  conven¬ 
tional  imaging  systems. 

1.  Ballistic  Imaging — Single  Point 

The  problem  of  determining  whether  an  object  lies 
along  the  line  of  sight  can  be  cast  as  a  statistical 
hypothesis  test  as  follows.  Given  a  received  photon 
count  X  at  the  sensor,  one  must  choose  between  two 
possible  situations.  The  first  situation  is  that  no  oc¬ 
cluding  object  exists  in  the  path  between  the  laser 
and  the  sensor  (H0).  The  alternative  is  that  there  is  an 
occluding  object  along  the  line  of  sight  between  the 
laser  and  the  detector  (H,).  Semitransparent  objects 
can  be  considered  to  be  significantly  more  scattering 
than  the  medium  [31],  and  therefore  the  detector  will 
collect  an  attenuated  number  of  ballistic  photons 
(compared  with  H0)  along  with  the  noise  photons. 

Following  the  notation  of  Subsection  3.B,  we  define 
the  number  of  noise  photons  that  will  arrive  at  the 
detector  as  99?(NXe),  where  X  ~  9‘Hx)  is  a  Poisson- 
distributed  random  variable  with  mean  x-  Then,  the 
hypothesis  test  is  given  by  the  null  hypothesis  de¬ 
fined  as  H0 :  X  ~  9HNXr  +  NX,.)  and  the  alternate 
hypothesis  (object  exists)  defined  as  Hx :  X  ~  9f'(NX„ 
+  NXS  ).  As  the  mean  of  the  Poisson  distribution 
grows,  the  probability  distribution  tends  to  a  Gauss¬ 
ian,  e.g.,  averaging  many  repeated  Poissonian  trials 
(i.e.,  N  large)  results  in  a  Gaussian-distributed  sta¬ 


tistic.  Using  the  Anscombe  transformation  [31],  we 
obtain  the  following  relationship: 

X~3>(x)^2jz7l~N(2jx,  1), 

where  ..V((,  a2)  represents  Gaussian  distribution  with 
mean  £  and  variance  a2.  Defining  a  new  variable 
representing  the  Anscombe-transformed  statistic  X' 
=  2\X  +  3/8  the  hypothesis  test  becomes 

H0  :  X' ~  N(2sNXe  +  NXS,  l), 

Hi  :  X'  ~N(2,NXe  +  NX,inc,  l).  (Al) 

The  decision  test  is  now  defined  as  (X'  $  jjj  y')> 
where  a  user-specified  false  alarm  rate  (PFA)  deter¬ 
mines  the  value  of  the  threshold  (y')  such  that  P(X' 

<  y  |  h0)  ^  pFa- 

2.  Ballistic  Imaging — K  Points 

The  problem  now  is  modified  to  describe  an  imaging 
scenario  of  scanning  the  FOV  at  a  fixed  square  array 
of  (ft  X  I\  i  points.  This  results  in  a  multiple  hypoth¬ 
esis  testing  problem  (K  tests),  where  for  large  K  it 
puts  a  lower  bound  on  the  SNR  of  the  observation.  To 
boost  the  SNR,  one  could  use  spatial  aggregation  by 
averaging  over  a  number  of  observation  points.  This 
modifies  the  problem  to  averaging  neighborhoods  of 
points  in  an  area  measuring  \W  X  W,  W  <  K,  effec¬ 
tively  reducing  the  spatial  resolution  of  the  detection 
map  (image).  By  decreasing  the  spatial  resolution, 
this  also  decreases  the  variance  at  each  point,  modi¬ 
fying  the  decision  test  to 

H0  :  X'~x[29NXe+NXs, 

Hi  :  X'  ~  TXe  +  NXSinc,  (A2) 

This  test  is  under  the  assumption  that  the  averag¬ 
ing  window  will  contain  either  no  occluder  points  or 
all  occluder  points.  In  reality,  the  averaging  filter  will 
result  in  an  observed  point,  X'  ~  .;\r(cpA[M'  |  H0]  + 
(1  —  9  )E\X'  |  HJ,  1/IF),  where  9  is  the  fraction  of  the 
window  containing  nonoccluders,  and  E  is  the  ex¬ 
pected  value  operator.  Our  goal  is  to  find  the  lower 
bound  on  the  value  of  W  that  will  guarantee  an  over¬ 
all  false  alarm  rate  of  less  than  P  FA  >  and  we  only 
consider  the  ideal  case  (all  occluders  or  nonoccluders) 
in  our  calculations  in  order  to  obtain  closed-form  so¬ 
lutions. 

The  Bonferroni  correction  approach  is  a  con¬ 
servative  method  of  controlling  the  false  alarm  rate 
for  a  detection  problem  under  multiple  independent 
and  identically  distributed  tests  [33].  The  correc¬ 
tion  adjusts  the  threshold  for  each  individual  test  in 
order  to  satisfy  a  lower  (per  test)  false  alarm  rate 
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value  ( Pfa/K)  such  that  each  of  the  fixed  number  K 
points  in  the  array  (and  W-point  averaging  filter) 
satisfies  [PCX’  <  y'  |  H0)  <  PFA/K\.  With  <b(x)  as  the 
cumulative  distribution  function  of  the  Af( 0,  1) 
density  at  the  point  x,  this  results  in  [y'  < 
(l/iW)®-\PFA/K)  +  2 ^NXe  +  WQ.  To  give  a  satis¬ 
factory  observation,  we  also  bound  the  miss  probabil¬ 
ity  for  detecting  a  ballistic  photon  by  the  same 
modified  value  ( PFA/K)  such  that  [PCX’  >  y'  |  H,) 
<  Pfa/K],  Using  the  miss  bounds,  we  determine  the 
lower  bound  on  the  necessary  averaging  window  size 
(W)  to  image  a  fixed  P-point  array  as 


- 1 

1  l—1 

o 

-<b  1 

-2 

[2  ,NXe  +  NXs- 

JNXe  +  NXs 

nc. 

The  minimum  width  of  the  occluding  object  {wb) 
that  can  be  reliably  resolved  for  a  given  parameter¬ 
ized  turbid  medium  can  now  be  derived.  Using  the 
lower  bound  for  W  found  in  Eq.  (A3),  we  can  solve 
for  the  lower  bound  on  the  width  using  wb  = 
v(FOV  X  W)/K. 

3.  Conventional  Imaging  Analysis 
In  the  conventional  imaging  regime,  there  is  no  time¬ 
gating  mechanism  and  all  the  photons  that  reach  the 
detector  over  a  long  acquisition  time  will  be  observed 
[acquisition  time  (d/c)  =  direct  line-of-sight  flight 
time].  Therefore,  a  large  number  of  photons  sent 
through  the  medium  will  be  collected  by  the  detector. 
A  problem  occurs  here,  too — although  the  SNR  is  high 
due  to  the  large  number  of  photons,  the  average  num¬ 
ber  of  scattering  events  on  each  photon  collected  will 


(g)  (h)  (i) 

Fig.  17.  Simulation  experiments  with  FOV  =  50  m  X  50  m,  dinc  =  0.3  m,  and  p,  =  12.5-1  m  C  (a)-(c)  Ballistic,  Bonferroni,  and 
conventional  observations  at  d  =  350  m,  respectively.  (dMf)  Ballistic,  Bonferroni,  and  conventional  observations  at  d  =  400  m,  respec¬ 
tively.  (g)— (i )  Ballistic,  Bonferroni,  and  conventional  observations  at  the  critical  distane  d  =  dcritical  =  417  m,  respectively. 
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also  be  high.  As  the  number  of  scattering  events  in¬ 
creases  for  a  photon,  the  spatial  resolution  of  the 
occluding  object  will  degrade.  The  lack  of  spatial  in¬ 
formation  results  in  a  blurred  observation.  Using 
random-walk  theory  [34,35],  it  is  possible  to  solve  for 
the  minimum  width  of  an  occluding  object  that  is  re¬ 
liably  resolved  using  a  conventional  imaging  system. 
The  width  [36]  is  found  using  the  photon  mean  time  of 
flight  (At),  which  can  be  numerically  computed  as  a 
function  of  the  parameters  of  the  medium  (p,s,  \x„). 
The  modified  minimum  full  width  at  half-maximum 
(FWHM)  is  equal  to  wconv  =  0.94((Ai)c/ ps)1/z. 

4.  Optimal  Resolution  Trade-Offs 

Ideally,  one  should  choose  the  imaging  system  (bal¬ 
listic  or  conventional)  that  reliably  resolves  the 
smallest  possible  object  \w  =  min whj\.  The  de¬ 
cision  test  using  the  minimum  resolvable  sizes  de¬ 
rived  above  becomes 


conventional 

Wconv  wb  =>  0.94 

ballistic 


/(Aty 

^^conventional 

FOV  x  W 

\  Fs 

\ 

ballistic 

K 

(A4) 


Using  the  lower  bound  of  W  from  Eq.  (A3),  one  can 
solve  for  the  critical  distance  ( d  =  d,.rltical),  the  maxi¬ 
mum  distance  at  which  ballistic  still  offers  superior 
resolution  relative  to  conventional  imaging. 

As  an  illustrative  example,  we  considered  a  ballis¬ 
tic  scanning  experiment  at  ( K  =  2562)  points,  imaging 
a  50  m  X  50  m  FOV.  The  occluding  objects  were  as¬ 
sumed  to  be  circular  of  diameter  1.0,  2.0,  4.0,  10.0, 
and  20.0  m  each  of  thickness  dinc  =  0.3  m  and  \xinc 
=  12.5  m-1.  We  used  false  alarm  rate  PFA  =  0.05  and 
considered  a  heavy  fog  turbid  medium  with  \x,  = 
12.5  1  m  '.  Using  the  analysis  from  above,  dcriticai  = 
417  m,  which  is  consistent  with  earlier  results  we 
showed  in  Section  3.  Figure  17  shows  the  effect  of 
distance  on  the  ballistic  resolution,  illustrating  the 
captured  ballistic  and  diffused  (conventional)  images 
at  distances  of  d  =  350,  400,  417  m. 
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RESOLUTION  BOUNDS  AND  RECONSTRUCTION  TECHNIQUES 
FOR  TIME-RESOLVED  TRANSILLUMINATION  IMAGING 

BRIAN  ERIKSSON  AND  ROBERT  NOWAK* 


Abstract. 

Recent  technological  advances  now  enable  time-gated  acquisition  of  photons  at  very  fast  rates. 
This  allows  one  to  separate  scattered  and  unscattered  photons  by  temporal  gating,  a  process  termed 
time-resolved  transillumination  (TRT)  imaging.  TRT  imaging  opens  the  door  to  a  new  type  of 
scanning  through  turbid  (scattering)  media  such  as  soft  tissue  and  fog/smoke,  and  many  poten¬ 
tial  applications  in  bioimaging  and  surveillance.  This  paper  investigates  the  performance  tradeoffs 
between  TRT  imaging  and  conventional  (un-gated)  imaging.  On  the  one  hand,  time-gated,  first- 
arrival  photons  are  unscattered  and  therefore  provide  very  high  spatial  resolution.  But,  very  few 
photons  arrive  at  the  detector  without  scattering,  effectively  resulting  in  a  very  low  SNR.  On  the 
other  hand,  conventional  (un-gated)  imaging  is  based  on  all  photons  (scattered  and  unscattered), 
resulting  in  lower  spatial  resolution,  but  higher  SNR  (due  to  the  large  number  of  photons).  This 
paper  investigates  these  tradeoffs  using  a  decision-theoretic  approach  to  ascertain  bounds  on  the 
minimum  resolvable  occluding  object  size  with  and  without  time-gated  photon  acquisition.  The 
theoretical  predictions  are  validated  through  a  realistic  simulation  of  tumors  in  breast  tissue.  The 
paper  then  proposes  a  novel  Maximum  Likelihood  approach  to  TRT  image  reconstruction.  Using 
a  novel  Expectation-Maximization  algorithm,  multiple  snapshot  observations  of  time-gated  photons 
are  used  to  reconstruct  the  image.  This  allows  us  to  obtain  the  ”  best-of-both- worlds”  by  combining 
the  spatial  resolution  of  the  first-arrival  photons,  with  the  higher  SNR  provided  by  scattered  photons. 

Key  words.  Transillumination  imaging,  Time-gated  imaging,  Photon- limited  imaging,  Image 
reconstruction,  High-resolution  imaging 

AMS  subject  classifications.  62C05,  94A08 


1.  Introduction.  Recent  technological  advances  now  enable  time-gated  acqui¬ 
sitions  of  photons  at  very  fast  rates,  fast  enough  to  separately  collect  unscattered 
(first-arrival)  and  scattered  (later-arrival)  photons  in  transillumination  imaging  [1], 
We  refer  to  this  technology  as  time-resolved  transillumination  (TRT)  imaging.  TRT 
imaging  opens  the  door  to  a  new  type  of  imaging  through  turbid  (scattering)  media 
(e.g.,  soft  tissue,  fog/smoke).  The  ability  to  separately  detect  unscattered,  so-called 
ballistic,  photons  can  enable  much  higher  spatial  resolution  imaging  than  possible 
using  conventional  (un-gated)  imaging  devices,  and  this  has  many  potential  applica¬ 
tions  in  bioimaging  and  surveillance.  However,  a  tradeoff  exists,  since  the  number  of 
ballistic  photons  decays  exponentially  as  the  thickness/depth  of  the  turbid  medium 
increases.  Therefore,  the  high  resolution  information  that  is  available  is  often  in  a 
very  photon-limited,  low  SNR  regime.  This  paper  explores  a  novel  approach  to  TRT 
image  reconstruction  and  a  decision-theoretic  approach  towards  deriving  bounds  on 
achievable  spatial  resolution. 

In  more  detail,  the  TRT  imaging  problem  involves  photons  traveling  through  a 
turbid  medium  from  a  source  (usually  a  laser),  through  a  scattering  medium,  and 
then  onto  a  detector  plane  as  depicted  in  Fig.  1.1.  Photons  traveling  through  a 
scattering  medium  can  be  roughly  classified  into  three  types:  ballistic,  snake,  and 
diffuse.  Ballistic  photons  experience  no  scattering  and  thus  travel  in  a  direct  line 
of  sight,  arriving  first  at  the  detector.  Because  of  the  lack  of  scattering,  ballistic 
photons  retain  their  spatial  information  and  arrive  at  the  detector  plane  at  the  same 
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planar  location  as  sent  from  the  light  source.  In  the  intermediary  regime,  ’’snake 
photons”  experience  slight  scattering  through  the  medium,  this  scattering  will  cause 
these  photons  to  arrive  later  than  the  ballistic  photons  and  likely  in  a  slightly  different 
planar  location  than  the  light  source  (resulting  in  a  small  loss  of  spatial  resolution). 
The  third  class,  diffuse  photons,  experience  large  amounts  of  scattering  and  arrive 
at  the  image  plane  having  lost  most  of  their  point  of  origin  information.  Due  to  the 
large  number  of  scattering  events  through  the  medium,  the  diffuse  photons  will  travel 
the  farthest  distance  to  the  detector  plane,  and  therefore  will  arrive  after  the  snake 
and  ballistic  photons.  While  the  inherent  spatial  information  decreases  in  order  of 
ballistic,  snake  and  diffuse  photons,  the  number  of  photons  (and  hence  inherent  SNR) 
increases  in  the  same  order.  So  we  are  faced  with  high  resolution,  low  SNR  data 
at  one  extreme  (ballistic),  and  low  resolution,  high  SNR  data  at  the  other  (diffuse). 
Furthermore,  the  diffusion  and  SNR  parameters,  which  characterize  the  underlying 
point  spread  function  (PSF),  are  not  known  precisely  in  practice. 


Fig.  1.1.  Example  of  Photons  Through  a  Scattering  Medium 


1.1.  Related  Works.  The  majority  of  the  previous  work  in  the  area  of  transillu¬ 
mination  imaging  has  been  focused  on  finding  the  resolution  limits  on  an  observation 
using  the  diffuse  photons  [1] ,  [2] ,  [3] ,  [4] ,  [5] ,  [6] .  Relatively  little  work  has  been  performed 
on  characterizing  the  resolution  limits  of  the  ballistic  observation  (most  notably  [7] 
and  [1],  neither  of  which  take  detection  error  rates  into  careful  consideration).  To 
the  best  of  our  knowledge,  this  is  the  first  attempt  to  quantify  the  tradeoff  between 
ballistic  and  diffuse  photons  in  terms  of  achievable  spatial  resolution.  Previous  work 
has  been  performed  on  iterative  least-squares  methods  to  reconstruct  images  based 
on  time-resolved  turbid  medium  observations  ([8],  [9]).  This  paper  is  the  first  attempt 
at  using  a  maximum  likelihood  approach  to  estimate  turbid  medium  parameters,  and 
also  the  first  attempt  at  developing  an  EM  algorithm  for  transillumination  image 
reconstruction. 

1.2.  Organization.  This  paper  is  divided  into  three  parts.  Section  2  uses  a 
decision-theoretic  approach  to  derive  resolution  bounds  for  turbid  environments.  In 
Section  3  we  introduce  an  estimation-theoretic  approach  to  TRT  image  reconstruction. 
A  statistical  model  for  TRT  imaging  is  purposed,  and  then  an  EM  algorithm  using 
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a  penalized  maximum  likelihood  approach  is  described.  Finally,  in  Section  4  we 
will  simulate  transillumination  imaging  through  breast  tissue  and  examine  the  EM 
algorithm  image  reconstruction  results. 

2.  Resolution  Bounds  for  Conventional  and  Time-Resolved  Imaging. 

The  resolution  of  an  transillumination  imaging  system  can  be  quantified  in  terms  of 
the  smallest  occluding  object  that  can  be  reliably  detected  and  located.  To  address 
this  issue,  we  consider  the  following  imaging  process.  Assume  that  the  source  and 
detector  pair  scan  the  area  of  the  detector  plane  (A)  at  points  on  an  n  x  n  grid,  to 
create  a  pixelated  cross-section  of  the  detector  plane  area  A.  We  assume  that  the 
region  between  the  source  and  detector  is  a  homogenous  scattering  medium,  possibly 
containing  occluding  objects  embedded  in  the  region.  Our  goal  is  to  determine  the 
smallest  occluding  object  that  can  be  reliably  detected.  The  size  of  an  occluder 
is  measured  by  its  cross-sectional  diameter  ( w°  in  Figure  2)  and  its  thickness  (d°  in 
Figure  2).  We  adopt  a  cylindrical  object  model  so  that  the  scattering  properties  along 
the  line-of-sight  between  the  source  and  detector  are  spatially  uniform  (a  spherical 
object  model  could  also  be  considered,  but  because  its  thickness  would  vary  depending 
on  where  the  line-of-sight  intersects  the  sphere,  its  analysis  would  be  slightly  more 
complicated) . 

For  this  section,  only  the  two  extremes  of  the  imaging  regime  will  be  considered, 
the  ballistic  regime  where  only  early  arriving,  unscattered  photons  are  detected  using 
high-speed  gating  mechanisms,  and  the  conventional  imaging  regime  where  all  the 
photons  arriving  at  the  detector  are  utilized  (i.e.  no  time-gating).  We  characterize 
the  minimum  size  occluder  that  can  be  reliable  resolved  in  these  two  regimes,  and 
show  that  there  is  a  distinct  threshold  on  the  total  image  acquisition  time,  below 
which  conventional  imaging  offers  better  spatial  resolution  an  above  which  ballistic 
imaging  is  superior.  This  threshold  is  a  function  of  the  parameters  of  the  scattering 
medium  and  occluding  material. 


Scattering  Medium 


Source 

Point 


Detector 


Point 


Fig.  2.1.  Turbid  Medium  Model  Example 


A  turbid  medium  may  be  considered  as  any  scattering  material,  but  for  example 
purposes  throughout  the  paper  we  will  assume  it  to  be  a  homogeneously  scattering 
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section  of  human  tissue.  The  scattering  material  is  parameterized  by  the  number 
of  photon  scatters  per  unit  length  (the  scattering  coefficient)  //"* ,  the  fraction  of 
absorptions  per  length  (the  absorbtion  coefficient)  gL™,  and  the  length  of  the  medium, 
dm .  Between  the  source  and  detector  in  the  medium  there  may  exist  occluding  objects, 
which  could  represent  tumors  in  the  bioimaging  context.  For  the  purposes  of  our 
resolution  bound  analysis,  we  assume  that  the  occluders  are  located  in  the  plane 
midway  between  the  source  detector,  since  this  is  the  plane  in  the  medium  where 
the  variance  of  the  scattering  will  be  at  a  maximum  [2], [3].  To  completely  define  the 
environment,  one  must  take  into  consideration  the  scattering  properties  of  occluders 
(scattering  and  absorption  coefficients  denoted  by  fi°,  n°,  respectively).  In  general,  the 
numbers,  locations  and  geometrical  properties  of  occluding  objects  may  be  arbitrary, 
but  for  the  purposes  of  determining  the  minimum  resolvable  occluder  size,  we  assume 
a  cylindrical  occluder  of  thickness  d°  and  width  w°. 

We  first  consider  the  spatial  resolution  of  conventional  (un-gated)  imaging,  which 
can  be  determined  via  classical  random-walk  diffusion  techniques  and  depends  only 
on  the  properties  of  the  scattering  medium,  and  not  on  the  particular  properties  of  the 
occluders.  Then  we  consider  the  spatial  resolution  achievable  from  ballistic  photon 
imaging.  Since  the  numbers  of  detected  ballistic  photons  depend  on  the  properties 
of  both  the  scattering  medium  and  the  occluders,  this  analysis  is  somewhat  more 
complex  than  that  of  conventional  imaging.  We  cast  the  ballistic  imaging  problem 
as  a  multiple  testing  problem,  where  our  goal  is  to  reliable  decide  whether  or  not  an 
occluder  lies  along  the  line-of-sight  between  the  source  and  detector  at  each  of  the 
n  x  n  scan  locations. 

2.1.  Conventional  Imaging  Resolution.  In  the  conventional  imaging,  there 
is  no  time-gating  mechanism  and  all  the  photons  that  reach  the  detector,  over  a 
relatively  long  integration  time,  are  recorded.  This  results  in  a  high  SNR,  but  since 
most  of  the  photons  are  scattered,  the  spatial  resolution  is  limited.  In  effect,  the 
resulting  image  is  blurred,  but  not  very  noisy.  Methods  for  determining  the  blurring 
caused  by  scattering  vary  in  complexity.  Simple  linear  techniques  found  in  [4], [10] 
have  been  shown  to  fit  experimental  data  well  for  short  integration  times.  Due  to  the 
absorption  of  photons  traveling  long  distances,  these  linear  techniques  do  not  hold  as 
the  integration  time  tends  to  infinity.  More  complicated  methods  using  Monte  Carlo 
techniques  ([6]),  or  diffusion  theory  calculations  ([1], [5])  are  found  to  be  more  accurate 
when  considering  long  integration  times.  Another  method  accurate  with  respect  to 
long  integration  times  uses  random  walk  theory  in  [2], [3],  and  due  to  the  accuracy  and 
simplicity  of  the  method  we  will  consider  this  approach.  The  blurring  effect  caused 
by  scattering  can  be  modeled  as  a  convolution  with  a  Gaussian  kernel.  The  standard 
deviation  of  the  Gaussian  kernel  is  found  using  the  photon  mean-time-of- flight  (At), 
which  can  be  numerically  solved  ([2],  Eqn.  16)  using  the  parameters  of  the  medium 
(/r™,/z™).  The  full- width  half-maximum  of  the  Gaussian  kernel  is 


w 


O 

conv 


(2.1) 


which  provides  a  reasonable  measure  of  the  spatial  resolution  of  a  conventional  imaging 
system  operating  in  transillumination  mode  through  the  given  medium. 

2.2.  Ballistic  Imaging  Resolution.  To  begin,  let  us  first  consider  the  statistics 
of  detected  ballistic  photons  along  a  single  line-of-sight.  The  problem  of  determining 
whether  or  not  an  occluder  lies  along  the  line-of-sight  can  be  cast  as  a  statistical 
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hypothesis  test,  as  follows.  Given  the  ballistic  photon  count  X  at  the  detector,  one 
must  choose  between  two  hypotheses.  The  first  is  that  no  occluding  object  exists  in 
the  path  between  the  source  and  detector,  which  we  denote  by  TIq.  The  alternative 
is  that  there  is  an  occluder  along  the  line-of-sight,  denoted  by  H\.  Occluders  are 
considered  to  be  significantly  more  scattering  than  the  surrounding  medium  (e.g., 
tumors  in  human  tissue  [11]),  and  therefore  the  detector  will  collect  an  attenuated 
number  of  ballistic  photons  under  TL[  compared  to  7io). 

The  probability  of  a  photon  passing  through  the  scattering  medium  with  no  scat¬ 
tering  events  (i.e.  the  probability  of  a  photon  of  being  a  ballistic  photon)  is  exponen¬ 
tially  distributed  with  rate  parameter  (/z™  +  /z™)  drn .  Given  the  source  intensity  A, 
equal  to  the  number  of  photons  sent  through  the  medium  by  the  per  unit  time,  and 
assuming  that  no  occluder  lies  along  the  path  between  source  and  detector,  one  can 
determine  the  intensity  of  the  ballistic  photon  arrival  process  at  the  detector: 

A0  =  Aexp(-GC  +  /C  )dm)  (2.2) 

When  an  occluder  does  lie  along  the  path,  the  intensity  of  the  ballistic  photon  arrival 
process  is  attenuated: 

Az  =  Aexp(-(/z™  +  *C  )(dm-d°)-(ix°s+^a)d0)  (2.3) 

In  addition  to  ballistic  photons  due  to  the  source,  stray  “noise”  photons  due  to  back¬ 
ground  light  sources,  arriving  during  the  time-gate  interval  of  acquisition,  will  also  be 
detected.  Let  A&  denote  the  average  intensity  of  this  arrival  process. 

We  can  now  express  the  probability  distributions  of  the  photon  count  at  the 
detector,  under  the  two  hypotheses. 


H0  :X  ~V((\b  +  \0)t) 


(2.4) 


Hi:X~V{{ A6  +  Ai)t)  ,  (2.5) 

where  t  is  the  length  of  the  time-gate  interval  in  seconds  and  'P(A)  denotes  the  Poisson 
distribution  whose  probability  mass  function  is  given  by  p(X  =  to)  =  e~AAm/m\.  As 
the  mean,  A,  of  the  Poisson  distribution  grows,  the  probability  distribution  tends  to  a 
Gaussian.  Averaging  repeated  trials  (i.e.,  averaging  of  multiple  laser  pulses)  therefore 
results  in  a  Gaussian  distributed  statistic.  Using  the  Anscombe  Transformation  [12], 
we  can  obtain  approximately  Gaussian  distributed  statistics.  If  X  ~  'P(A),  then  the 

Anscombe  transformed  variable  2  X  +  |  is  asymptotically  (as  A  — >  oo)  Gaussian 
distributed  as  A f  ^2-y/A,  1^  (where  Af  (/i,  a2)  is  a  Gaussian  distribution  with  mean  /./ 

and  variance  a2).  In  fact,  the  approximation  to  Gaussian  is  quite  accurate  even  when 
A  is  relatively  small  (e.g.,  A  =  20).  Thus,  we  transform  the  photon  count  X  to  obtain 

a  new  statistic  X'  =  2y ' X  +  |.  Letting  T  =  kt  denote  the  total  acquisition  time, 
over  k  trials  with  t  units  of  time  per  trial,  the  resulting  Gaussian  hypotheses  are  given 
by: 


na:X'~M  (2^(Ab  +  A0 )T,  l) 
Hi-.X'-M  ^2 a/ (Ai  +  Ai)  T,  l) 


(2.6) 

(2.7) 
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The  decision  test  is  now  defined  as  55^*  y'j  .  The  probability  of  error  of  this  test 
can  be  expressed  in  terms  of  the  Gaussian  error  function  <f>(f)  :=  J)°°  -^e-4  ^2>  as 

P  =  2$(v/(Ab  +  A0)T-  (Ab  +  Ai)T^  =  2$  (%/TAa) 


where  A  a  —  a/A;,  +  Ao  —  \/Ab  +  Ax . 

Now  recall  that  our  goal  is  to  create  a  pixelated  image  atnxn  points  over  the  FOV. 
In  other  words,  the  testing  procedure  above  will  be  repeated  at  n2  different  locations. 
Whenever  a  statistical  test  is  repeated  multiple  times,  the  chance  that  one  or  more  is 
in  error  increases  with  the  number  of  trials.  Our  goal  is  to  guarantee  a  low  probability 
of  error  at  all  n2  locations  (i.e.,  to  control  the  probability  that  an  erroneous  decision 
is  made  at  one  or  more  of  the  pixel  locations).  This  probability  is  bounded  by  the 
sum  of  the  probability  of  error  at  each  location.  That  is,  the  overall  probability  that 
one  or  more  pixels  is  in  error  is  bounded  by  n2p.  This  bound  is  called  the  Bonferroni 
or  union  bound.  Thus,  if  we  desire  an  overall  probability  of  error  pe,  then  we  require 
that 


n  <  VPe/P 


1/2 


The  minimum  width  of  an  occluder  that  can  be  reliably  resolved  for  a  given  turbid 
medium  can  now  be  derived.  Using  the  upper  bound  for  n  above,  we  can  solve  for  the 

lower  bound  on  the  width  using  w^aUistic  =  \j~^i  (where  A  is  the  area  of  the  detector 
plane) : 


^ ballistic  \/ '  AJfl 


> 


1/2 


(2.8) 

(2.9) 


This  equation  shows  how  the  minimum  resolvable  occluder  size  depends  on  the  pa¬ 
rameters  of  the  imaging  system. 

This  resolution  can  be  contrasted  with  that  of  the  conventional  imaging  system 
given  by  (2.1).  The  critical  acquisition  time,  Tcritical  (per  pixel),  at  which  the  resolution 
of  the  ballistic  system  w%aUistic  is  equal  to  that  of  the  conventional  system  w°onv ,  is 
given  by 


TC1 


A  a  <f)_1 


Pe  (At)  C 
2  AW 


2 


(2.10) 


At  acquisition  rates  below  Tcritical,  conventional  imaging  provides  better  resolution 
than  ballistic  imaging,  and  the  converse  is  true  at  faster  rates  above  Toritical.  Deter¬ 
mination  of  Tcritical  requires  specification  of  the  parameters  of  the  turbid  medium  and 
occluding  material,  and  the  allowable  error  level  in  the  ballistic  case  (e.g.,  pe  =  0.05). 
We  will  examine  the  critical  acquisition  time  in  examples  in  the  following  sections. 

3.  TRT  Image  Reconstruction.  The  analysis  above  demonstrates  that  bal¬ 
listic  imaging  can  provide  statistically  reliable,  higher  resolution  images  than  conven¬ 
tional  imaging  in  certain  situations.  In  this  section,  we  develop  a  novel  approach  to 
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TRT  image  reconstruction  that  takes  advantage  of  both  the  high  resolution  of  bal¬ 
listic  photons,  and  the  higher  SNR  associated  with  conventional  imaging,  leveraging 
the  best  of  both  regimes. 

The  TRT  image  reconstruction  problem  can  be  posed  as  a  statistical  inverse 
problem  (a  particular  form  of  photon-limited  image  reconstruction),  but  to  the  best  of 
our  knowledge  our  work  here  is  the  first  to  formulate  it  as  such.  Due  to  uncertainties  in 
the  blurring  kernel  (dependent  on  usually  unknown  parameters  of  the  turbid  medium) , 
the  TRT  problem  is  somewhat  related  to  blind  image  deconvolution  (BID).  While 
BID  is  a  well-studied  problem,  there  are  several  unique  aspects  in  the  TRT  imaging 
problem  that  make  it  quite  different  from  standard  BID  problems.  In  particular, 
the  distinctive  features  of  TRT  imaging  include  the  photon-limited  nature  of  the 
data,  the  time-gated  data  acquisition  which  in  effect  yields  information  at  multiple 
spatial  resolutions  that  can  be  fairly  well  characterized  via  a  diffusion  equation,  and 
most  importantly  the  availability  of  “unblurred”  data  corresponding  to  the  ballistic 
photons.  For  these  reasons  “off-the-shelf”  BID  algorithms  are  not  directly  applicable 
to  TRT  image  reconstruction. 

3.1.  A  Statistical  Model  of  TRT  Imaging.  The  basic  statistical  model  we 
propose  for  TRT  imaging  through  homogeneous  turbid  media  is  as  follows.  Assume 
that  we  have  k  time-resolved  “snapshots”,  each  acquired  over  disjoint  time  intervals 
. . . ,  T^k\  with  denoting  the  ballistic  time  interval.  Assume  that  these  inter¬ 
vals  form  a  uniform  partition  of  the  overall  observation  interval,  denoted  T  >  0.  Let 
X^\ . . . ,  denote  the  photon  data  acquired  in  each  interval,  respectively.  Specif¬ 
ically,  the  data  X M  are  acquired  in  the  form  of  an  n  x  n  pixel  image,  and  for  our 
mathematical  exposition  we  assume  that  the  columns  of  this  image  are  “stacked”  to 
form  an  n2  x  1  vector.  Each  entry  in  is  simply  the  number  of  photons  detected 
at  the  corresponding  pixel  location  during  the  time  interval  ;(-T) .  We 

assume  that  each  image  is  Poisson  distributed  according  to  the  following  model: 

X®  ~  Poisson(awP(i)A),  i  =  1, . . . ,  k,  (3.1) 

where  A  denotes  the  underlying  n2xl  image  intensity  function  we  wish  to  reconstruct, 
P M  denotes  the  n2  x  n2  photon  transition  matrix  (from  the  emission  (source)  plane 
to  the  detector  plane)  reflecting  the  diffusion  process  governing  photon  migration 
during  each  time  interval,  and  >  0  is  a  scalar  gain  factor  reflecting  the  varying 
abundances  of  detected  photons  in  the  differing  regimes,  from  ballistic  to  diffuse. 

The  transition  matrices  are  functions  of  time  and  a  scalar  diffusion  bandwidth 
parameter  denoted  by  er2 .  Using  the  elementary  linear  models  of  the  photon  migration 
and  diffusion  process  found  in  [4, 10] ,  we  model  PW  as  a  linear  shift-invariant  Gaussian 
blurring  kernel  with  variance  proportional  to  a2t^\  where  denotes  the  midpoint 
of  the  i-tli  acquisition  time-interval.  Thus,  the  transition  matrices  are  parametric 
functions  of  the  form  pM  =  P(cr2t W),  where  only  the  variance  coefficient  er2  needs 
to  be  estimated.  Also,  due  to  the  linear  shift-invariant  nature  of  the  model,  the 
resulting  blurring  process  is  a  convolution,  and  thus  the  forward  solution  can  be 
rapidly  computed  via  the  Fast  Fourier  Transform  (FFT)  algorithm  (a  property  we 
will  exploit  in  our  reconstruction  algorithms). 

We  assume  that  A  is  normalized  such  that  the  intensities  in  A  sum  to  one  (i.e., 
]P'l=1  A (j)  =  1).  Also  the  transition  matrices  are  normalized  so  that  //W  =  p also 
satisfies  Y^j=i  ^  U)  =  1-  With  these  normalizations,  it  is  easy  to  verify  that  the 
total  intensity  (expected  total  number  of  detected  photons)  of  the  snapshot  X^  is  the 


gain  Q!j,  i  =  1, . . . ,  fc.  Furthermore,  we  adopt  the  convention  that  the  ballistic  image 
resolution  is  the  finest  spatial  resolution  available  and  assume  that  P W  is  the  n2  x  n2 
identity  matrix.  We  also  assume  that  the  images  acquired  in  each  time-interval  are 
statistically  independent,  and  so  the  joint  distribution  function  of  the  entire  kn2  x  1 
data  record,  obtained  by  stacking  the  snapshots  into  a  single  column  vector 

-  xd)  ' 
x=  i 

xw 


is 


X  ~  Poisson(PA),  (3.2) 

where  P  is  the  kn 2  x  n 2  transition  matrix  obtained  by  stacking  the  matrices 
aWpW,...,aWp(fe);  i.e., 


P 


a«P C1)  ' 

a{k)  p(k) 


Let  us  contrast  the  above  model  of  the  TRT  imaging  system  with  a  conventional 
imaging  system  in  which  the  photon  detections  are  not  time-resolved.  In  this  case,  we 
acquire  the  aggregated  photon  image  Xa  =  Af^  +  •  •  •  +  X^k\  which  obeys  the  model 

Xa  ~  Poisson(PaA),  (3.3) 

where  Pa  =  aWp(*)..  We  will  see  that  the  extra  information  available  in  the 

time-resolved  photon  acquisition  can  significantly  improve  our  ability  to  estimate  the 
underlying  intensity  A. 

3.2.  Maximum  Likelihood  Reconstructions.  Both  the  TRT  imaging  model, 
X  ~  Poisson(PA),  and  the  conventional  imaging  model,  Xa  ~  Poisson(PaA),  are  in 
the  form  of  the  basic  statistical  model  usually  considered  in  photon-limited  imaging 
problems  [13].  If  the  parameters  of  the  imaging  system  (i.e.,  the  gains  a,  and  the 
diffusion  parameter  a2)  are  known,  then  standard  techniques  can  be  applied  to  find 
the  Maximum  Likelihood  Estimate  (MLE)  of  A. 

To  define  the  MLE,  let  p(Af|A)  denote  the  likelihood  function  of  the  data  (i.e., 
the  Poisson  probability  mass  function  evaluated  at  X  and  viewed  as  a  function  of  A) . 
The  MLE  is  the  solution  to  the  following  optimization: 

Amle  =  argmaxp(X|A)  (3.4) 

More  specifically,  taking  the  logarithm  of  the  likelihood  function  (a  monotonic  trans¬ 
formation),  we  have 


k 

Xmle  =  argmax  ^  (a(i)PwA  +  X®  loga(i)P(i)A  -  log (3.5) 

i—1 

There  is  no  closed-form  solution  for  the  Xmle ,  but  the  log-likelihood  function  is  con¬ 
vex  in  A.  Thus,  numerical  optimization  techniques  can  be  readily  applied.  One  very 
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popular  and  effective  optimization  procedure  is  the  so-called  Expectation-Maximization 
(EM)  algorithm  [13].  The  EM  algorithm  is  an  iterative  procedure  that  is  guaranteed 
to  produce  a  sequence  of  iterates,  A^1),  A^2\  . . . ,  that  have  non-decreasing  likelihood 
values  (i.e.,  each  iteration  of  the  algorithm  does  not  decrease  the  likelihood  value). 
The  EM  algorithm  can  work  very  well  in  practice,  but  usually  requires  an  ad  hoc 
stopping  criterion,  since  the  MLE  solution  can  very  noisy  in  severely  photon-limited 
and/or  ill-posed  cases.  Stopping  the  algorithm  prior  to  convergence  provides  an  indi¬ 
rect  smoothing  of  the  solution. 

A  more  desirable  alternative  to  the  early  stopping  of  the  EM  algorithm  is  to 
tackle  the  undesirable  aspects  of  the  MLE  solution  more  directly.  This  is  usually 
accomplished  by  augmenting  the  likelihood  criterion  with  an  additive  penalty  func¬ 
tion  which  discourages  solutions  that  are  rough  and  irregular.  The  generic  form  of 
penalized  solutions  is 

A  =  argmax  [log(p(X|A)  +  pen(A)]  (3.6) 

where  pen(A)  >  0  is  a  penalty  function  that  assigns  small  values  to  smooth  A  and 
larger  values  to  non-smooth  A.  The  notion  of  smoothness  depends  on  the  nature  of 
the  penalty,  and  many  possibilities  have  been  studied  in  the  literature  (see  [14]  and 
the  references  therein).  In  particular,  we  advocate  the  use  of  multiscale  complexity 
penalty  [15,  14],  which  tends  to  favor  solutions  that  are  piecewise  smooth.  The  mul¬ 
tiscale  complexity  penalty  leads  to  a  denoising  process  akin  to  wavelet  thresholding 
methods,  but  particularly  adapted  to  Poisson  data  [16,  15,  17,  14].  Figure  3.1  illus¬ 
trates  the  denoising  capability  of  the  multiscale  complexity  penalty.  In  this  example, 
a  single  photon-limited  snapshot  was  acquired,  without  blur  (i.e.,  the  corresponding 
photon  transition  matrix  is  the  identity  operator).  The  denoising  process  is  carried 
out  via  a  computationally  efficient  multiscale  tree-pruning  process  (akin  to  wavelet 
coefficient  thresholding)  that  requires  0(n2)  operations  for  an  n  x  n  image. 


(a)  (b)  (c) 

Fig.  3.1.  Night  sky  image  (a)  Original  Image,  (b)  Poisson  photon  count  image,  (c)  MPLE 
Denoised  Image,  resulting  from  the  MPLE  proposed  in 


The  multiscale  complexity  penalty  can  also  be  easily  integrated  into  an  EM  al¬ 
gorithm,  by  replacing  the  conventional  (MLE)  M-step  with  a  denoising  step  (see  [17] 
for  further  details).  Iterating  this  modified  EM  algorithm  to  convergence  produces 
a  Maximum  Penalized  Likelihood  Estimate  (MPLE).  To  the  best  of  our  knowledge, 
the  quality  of  the  MPLE  solution,  in  terms  of  squared  error,  is  state-of-the-art  [17], 
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and  so  we  will  exclusively  use  this  reconstruction  process  here.  Each  iteration  of  the 
EM  algorithm  follows  a  two-step  procedure.  In  the  E-step,  the  conditional  expecta¬ 
tion  of  the  so-called  “complete  data”  is  computed.  The  complete  data  is  an  n2  x  n2 
array  in  which  the  i.  j  element  is  equal  to  the  number  of  photons  emitted  from  posi¬ 
tion  i  and  detected  at  position  j.  Of  course,  we  only  observe  the  detected  photons, 
without  knowledge  of  where  they  originated  from,  and  so  the  complete  data  must 
be  estimated  based  on  the  photon  detections  X  and  the  current  iterate  of  A.  The 
E-step  requires  computing  matrix  vector  products  with  the  transition  matrices  PM 
(see  pseudocode  below).  Since  the  transition  matrix  operation  is  a  linear  convolu¬ 
tion,  these  products  can  be  rapidly  computed  via  the  Fast  Fourier  Transform  (FFT) 
algorithm.  From  the  complete  data,  we  can  compute  the  photon  emission  counts  by 
marginalization  (providing  data  as  if  no  diffusion /blurring  process  occurred).  The 
exact  (not  estimated)  photon  emission  counts  gives  the  MLE  for  A.  Thus,  the  stan¬ 
dard  M-step  (corresponding  to  the  MLE)  simply  computing  the  estimated  photon 
emission  counts  by  marginalizing  the  expectation  of  the  complete  data.  The  modified 
M-step  (corresponding  to  the  MPLE)  is  obtained  by  applying  the  multiscale  complex¬ 
ity  penalty  denoising  process  [15,  14]  to  the  estimated  photon  emission  counts.  The 
pseudocode  below  outlines  the  steps  of  the  EM  algorithm.  The  algorithm  produces  a 
sequence  of  iterates,  and  the  corresponding  sequence  of  penalized  log-likelihood  values 
is  monotonic,  non-decreasing  (see  [17]  for  further  details). 


Algorithm  1  -  MPLE  EM  Algorithm  with  Known  Parameters 


Initialize:  Set  iteration  counter  j  =  0  and  A (°\m)  =  1  /n2,  m=l,. . .  ,n2. 
E-Step:  Compute  the  conditional  expectation  of  the  complete  data  (c.f.  [13]): 

*«>(<,  m)  =  Tx>->(e)  f  W°(‘)pll|«,m) 


M-Step: 

1.  Compute  the  expected  photon  emission  counts  via  marginalization: 

n2 

(to)  =  ^ (£,m)  ,  i  =  l,..,,k 
e=i 

2.  Apply  the  multiscale  complexity  penalized  denoising  procedure  to 

to  obtain  updated  iterate  A^+1h  The  Poisson  denoising  procedure  was 
first  proposed  in  [15],  and  the  implementation  employed  in  this  work  is 
a  translation  invariant  version  developed  by  R.  Willett.  Code  is  on-line  at  [18]. 

Repeat:  Set  j  =  j  +  1.  Repeat  E-  and  M-steps  until  convergence. 


3.3.  Adapting  to  Unknown  Turbid  Media.  One  of  the  major  challenges 
in  practical  imaging  problems  is  that  the  characteristics  of  the  turbid  medium  are 
usually  not  known  precisely.  In  particular,  the  values  of  the  parameters  {a/®-*}  and  a2 
are  unknown  and  therefore  must  be  estimated  along  with  A.  We  can  formulate  this 
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as  a  joint  MLE  problem,  seeking  to  find  values  of  the  parameters  and  A  which  jointly 
maximize  the  Poisson  likelihood  function  (or  penalized  likelihood).  At  first  glance,  it 
may  appear  that  this  joint  MLE  problem  may  be  intractable,  but  it  turns  out  to  have 
a  rather  simple  solution. 

The  basic  solution  approach  is  as  follows.  First,  the  MLEs  of  the  gain  factors 
can  be  computed  separately  from  A  and  cr2  due  to  the  following  observation. 
Consider  the  statistics 


e=i 

These  statistics  are  simply  the  total  photon  counts  for  each  image.  Due  to  the  nor¬ 
malization  of  A  and  the  matrices  {PW}  in  our  model,  it  follows  that 

gW  ~  Poisson(a^),  i  =  l,...,k. 

It  is  well-known  that  the  conditional  distribution  of  given  S W  is  multinomial  with 
parameter  //W  =  P^A  (see  [15]).  Therefore,  the  likelihood  factorizes  into  Poisson 
factors,  each  involving  one  pair  (aW,SW),  and  multinomial  factors,  each  involving  A 
and  one  triple  (PW,XM,  S^).  Consequently,  the  MLEs  of  the  gain  factors  {a^}  can 
be  obtained  separately  from  A  and  cr2  and  are  given  by  the  simple  formula  aW  =  , 

i  =  1, . . . ,  k.  Now  recall  that  the  matrices  {P^}  are  parametric  in  a2.  To  find  the 
MLEs  of  A  and  a2  we  consider  a  range  of  candidate  values  for  cr2  and  for  each  one  we 
use  the  EM  algorithm  described  above  (with  each  a W  set  to  its  MLE  aW)  to  compute 
the  MPLE,  denoted  by  A(er2),  and  the  corresponding  maximum  penalized  likelihood 
value,  denoted  L(a2).  This  can  be  done  exhaustively  over  a  discretized  range  of  a2 
values,  by  performing  the  EM  algorithm  for  each  discretized  value  of  a2 .  The  value  of 
cr2  that  results  in  the  highest  penalized  likelihood  value  L(a2)  is  considered  the  MLE 
of  tr2,  denoted  by  a2.  The  MPLE  of  A  is  then  A(cr2).  Another,  less  computationally 
expensive  option,  is  to  incrementally  adjust  the  value  of  cr2  at  each  iteration  of  the 
EM  algorithm  (see  pseudocode  below  for  details) .  This  iterative  method  of  estimation 
is  performed  in  the  experimental  portion  of  this  paper  and  was  repeatedly  found  to 
converge  to  the  true  value  of  cr2. 

4.  Experimental  Results.  We  now  present  a  simulation  study  of  imaging  ma¬ 
lignant  breast  tissue  through  healthy  breast  tissue.  For  this  simulation,  a  medium  with 
a  100mm  x  100mm  detector  plane  area  (A)  and  length  dm  =  10mm  is  defined  with 
circular  occluding  objects  of  diameter  2mm,  4mm,  8mm,  20mm,  and  40mm.  From 
[11],  we  use  the  scattering  and  absorption  coefficients  of  the  two  tissue  types  (healthy 
breast  tissue,  malignant  tissue).  Using  a  specified  thickness  of  the  occluding  tumor 
(<P=2.5mm)  and  false  alarm  rate  (pe  =  0.05),  we  use  Equation  2.10  and  find  Tcritical  to 
equal  154ms.  We  simulate  a  series  of  ten  snapshot  observations  through  the  medium 
using  four  different  (per  snapshot)  acquisition  time  values  (T  =1,10, 100, 154ms). 

The  first  of  the  ten  snapshot  images  is  the  ballistic  observation  (as  seen  in  Fig¬ 
ure  4.1(A)  for  T=lms,  and  Figure  4.2(A)  for  T  =  Tcritical= 154ms,  generated  using  the 
formulas  in  Equations  2. 4, 2. 5),  and  the  last  is  a  heavily  blurred  diffuse  observation 
(Figure  4.1(B)  and  Figure  4.2(B)  generated  using  Equation  2.1).  A  single-snapshot 
version  of  the  MPLE  EM  algorithm  (i.e.,  with  k  =  1)  is  used  to  reconstruct  the 
ballistic  and  diffuse  observation  images  (Figure  4.1(C),  4.1(D)  for  T=lms,  and  Fig¬ 
ure  4.2(C),  4.2(D)  for  T  =  Tcritical= 154ms). 
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Algorithm  2  -  MPLE  EM  Algorithm  with  Unknown  Parameters 


Initialize: 

1.  Set  iteration  counter  j  =  0, 

2.  =  1/n2,  m=l,. . .  ,n2, 

3.  A«(f),  i  =  1, . . .  ,k. 

4.  Set  (tr^)2  =  0.164  (^),  an  elementary  estimate  of  the  sigma  parameter 
based  on  system  characteristics  (from  [10]) 

5.  Set  Act  =  i(cr(°))2, 


E-Step:  Set  P W  =  Compute  the  conditional  expectation  of  the  com¬ 

plete  data  (c.f.  [13]): 


Z U)  (£,  to) 


J2x(i)  w 

i= 1 


A ^  (to)  aWpM  (£,m) 
EtiAW  (m)a«pM  (*,m) 


M-Step: 

1.  Compute  the  expected  photon  emission  counts  via  marginalization: 

n2 

(to)  =  ^  (£,m)  ,  «  =  1, . . . ,  k 

e- i 

2.  Apply  the  multiscale  complexity  penalized  denoising  procedure  [18]  to 
to  obtain  updated  iterate  A^+1K 

3.  Set  (<j(j+1))2  =  argmax  L(cr2). 

CT2  =  {(ct0))2_A,t,((t0))2,(CT0))2+A(t} 

4.  Set  Aff  :  a2«. 


Repeat:  Set  j  =  j  +  1.  Repeat  E-  and  M-steps  until  convergence  and  record 


The  multiple  snapshot  EM  algorithm  (from  Section  3  is  used  to  fuse  the  data  from 
all  ten  snapshots.  Figure  4.1(E)  and  Figure  4.2(E)  show  the  image  reconstruction 
based  on  the  multiple  snapshot  data  given  exact  knowledge  of  all  attenuation  factors 
and  blurring  variances.  Using  the  adaptive  MLE  scheme  described  in  Section  3- 
3.3,  the  attenuation  factors  and  blurring  variance  were  jointly  estimated  along  with 
original  image  A,  with  the  resulting  image  reconstruction  shown  in  Figure  4.1(F)  and 
Figure  4.2(F).  Table  4.1  summarizes  the  reconstruction  errors  as  compared  with  the 
original  image  A.  As  expected,  the  combined  oracle  reconstruction  performs  best  (in 
terms  of  pSNR),  with  the  combined  reconstruction  with  ML  estimation  of  the  gain 
and  diffusion  parameters  close  behind. 

The  strength  of  the  EM  algorithm  can  be  observed,  as  the  EM  reconstructed  image 
will  tend  to  the  ballistic  reconstruction  at  T  >  Tcritical,  while  the  EM  reconstructed 


(A)  Ballistic  Obs. 


(B)  Diffuse  Obs. 


(C)  Ballistic  Reconst. 


(D)  Diffuse  Reconst. 


Fig.  4.1. 

$  =2.5  mm,  T  = 


(E)  Combined  Reconst. 


(F)  Unknown  Reconst. 


EM  Reconstruction  with  50  iterations,  10  images,  FOV  =  100mm  x  100mm  x  10mm, 
1  ms 


T 

Ballistic  Reconst. 

Diffuse  Reconst. 

Oracle 

Unknown 

1ms 

34.26 

39.33 

40.93 

40.43 

10ms 

38.09 

39.21 

41.08 

40.60 

100ms 

41.79 

39.26 

44.17 

43.85 

154ms 

42.31 

39.35 

44.88 

44.66 

Table  4.1 

Summary  of  TRT  image  reconstruction  Peak  Signal-to- Noise  Ratio  (pSNR  in  dB). 


image  will  tend  to  the  diffuse  reconstruction  at  T  <  Tcritical.  This  follows  the  analysis 
in  Section  2.  The  conservative  nature  of  the  Bonferroni  Correction  analysis  can  be 
seen  in  Table  4.1,  as  the  ballistic  observation  reconstruction  does  outperform  (in  a 
pSNR  sense)  the  diffuse  observation  reconstruction  for  acquisition  times  T  shorter 
than  the  Tcritk.al,  although  the  Tcritical  value  does  appear  to  be  of  the  correct  order  of 
magnitude. 

5.  Conclusions.  This  paper  presented  the  resolution  bounds  and  proposed  a 
novel  MLE  reconstruction  algorithm  for  TRT  imaging.  Using  the  decision-theoretic 
analysis  approach,  it  is  shown  that  for  a  given  turbid  medium,  the  smallest  reliably 
resolved  object  (in  terms  of  the  width  w°)  for  both  the  conventional  imaging  and 
ballistic  regimes,  can  be  derived.  The  image  reconstruction  algorithm  is  based  on 
a  combination  of  the  EM  algorithm  and  Poisson  denoising.  Our  simulation  study 
demonstrates  the  potential  of  our  approach,  in  particular  indicating  the  added  benefit 
of  optimally  fusing  ballistic  and  diffuse  photon  data  using  the  proposed  EM  algorithm. 
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(A)  Ballistic  Obs. 


(D)  Diffuse  Reconst. 


(B)  Diffuse  Obs. 


(E)  Combined  Reconst. 


(C)  Ballistic  Reconst. 


(F)  Unknown  Reconst. 


Fig.  4.2.  EM  Reconstruction  with  50  iterations,  10  images,  FOV  =  100mm  x  100mm  x  10mm, 
df  =  2.5  mm,  T  =  T critiCai  =  154 ms 
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ABSTRACT 

We  introduce  a  class  of  robust  non-parametric  estimation  meth¬ 
ods  which  are  ideally  suited  for  the  reconstruction  of  signals  and 
images  from  noise-corrupted  or  sparsely  collected  samples.  The 
filters  derived  from  this  class  are  locally  adapted  kernels  which 
take  into  account  both  the  local  density  of  the  available  samples, 
and  the  actual  values  of  these  samples.  As  such,  they  are  automat¬ 
ically  steered  and  adapted  to  both  the  given  sampling  “geometry”, 
and  the  samples’  “radiometry”.  As  the  framework  we  proposed 
does  not  rely  upon  specific  assumptions  about  noise  or  sampling 
distributions,  it  is  applicable  to  a  wide  class  of  problems  includ¬ 
ing  efficient  image  upscaling,  high  quality  reconstruction  of  an 
image  from  as  little  as  15%  of  its  (irregularly  sampled)  pixels, 
super-resolution  from  noisy  and  under-determined  data  sets,  state 
of  the  art  denoising  of  images  corrupted  by  Gaussian  and  other 
noise,  effective  removal  of  compression  artifacts;  and  more. 

Index  Terms —  Inverse  problem,  image  reconstruction,  piecewise 
polynomial  approximation,  nonlinear  estimation 

1.  INTRODUCTION 

Image  processing  methods  have  been  exploited  through  the  years 
to  improve  the  quality  of  digital  images.  Many  of  the  popular  im¬ 
age  processing  tools  have  a  limited  scope  of  use;  some  can  only 
be  employed  as  denosing  methods,  while  application  of  others 
are  limited  to  upscaling  regularly  sampled  data.  Moreover,  such 
methods  estimate  the  underlying  signal  based  on  certain  assump¬ 
tions  on  data  and  noise  models,  a  common  example  of  which  is 
modeling  the  noise  as  pure  additive  i.i.d.  Gaussian.  Although 
such  limiting  assumptions  facilitate  the  design  of  optimal  meth¬ 
ods  for  a  certain  type  of  data,  in  real  situations  when  the  data  and 
noise  models  do  not  faithfully  describe  the  measured  signal,  the 
performance  of  such  non-robust  methods  significantly  degrades 
[1]. 

Classical  parametric  image  processing  methods  rely  on  a  spe¬ 
cific  model  of  the  signal  of  interest,  and  seek  to  compute  the  pa¬ 
rameters  of  this  model  in  the  presence  of  noise.  In  contrast  to 
the  parametric  methods,  non-parametric  methods  rely  on  the  data 
itself  to  dictate  the  structure  of  the  model,  in  which  case  this  im¬ 
plicit  model  is  referred  to  as  a  regression  function  [2].  We  pro¬ 
mote  the  use  and  improve  upon  a  class  of  non-parametric  meth¬ 
ods  called  kernel  regression  [3],  which  generalizes  some  recently 
presented  methods  namely,  normalized  convolution  [4],  bilateral 
filter  [5,  6],  and  moving  least-squares  [7]. 

This  work  was  supported  in  part  by  DARPA/AFOSR  Grant  FA9550-06-1- 
0047;  by  AFOSR  Grant  F49620-03- 1-0387,  and  by  the  National  Science  Founda¬ 
tion  Science  and  Technology  Center  for  Adaptive  Optics,  managed  by  the  Univer¬ 
sity  of  California  at  Santa  Cruz  under  Cooperative  Agreement  No.  AST-9876783. 


The  main  advantage  of  the  presented  regression  method  is 
that  it  is  a  generic  framework  enabling  direct  use  in  a  variety  of 
applications,  from  single  frame  denoising  to  multi  frame  super¬ 
resolution  [3].  Moreover,  this  method  produces  better  and  more 
stable  results  comparing  to  the  state  of  the  art  methods  in  the 
literature,  as  it  is  robust  to  modeling  errors  and  data  outliers. 

This  paper  is  organized  as  follows.  Section  2  is  a  brief  intro¬ 
duction  to  the  notion  of  adaptive  kernel  regression  and  the  novel 
concept  of  using  weighted  l\  norm  penalty  term  in  the  kernel 
regression  framework.  Section  3  extends  and  generalizes  the  pre¬ 
vious  related  methods  to  derive  the  details  of  the  proposed  robust 
regression  method,  focusing  on  appropriate  choices  for  the  ker¬ 
nel  function.  Simulation  results  are  presented  in  Section  4,  and 
Section  5  concludes  this  paper. 

2.  DATA -ADAPTED  KERNEL  REGRESSION 

We  treat  the  2-D  estimation  problem  where  the  measured  data  yi 
at  the  position  X2 i]T  is  given  by 

yt  =  z(xi)  +  Ei,  i  =  1,  •••  ,P,  (1) 

where  z(-)  is  the  (hitherto  unspecified)  regression  function  (i.e. 
an  unknown  image)  to  be  estimated,  P  is  the  number  of  measured 
pixels,  and  efs  are  the  independent  and  identically  distributed 
noise  values  (with  otherwise  no  particular  statistical  distribution 
assumed). 

While  the  specific  form  of  z(-)  may  remain  unspecified,  if  we 
assume  that  it  is  locally  smooth  to  some  order  N,  then  in  order 
to  estimate  the  value  of  the  function  at  any  given  point  x,  we  can 
rely  on  a  generic  local  expansion  of  the  function  about  this  point. 
Specifically,  if  x  is  near  the  sample  at  x^,  we  have  the  IV-term 
Taylor  series 

z(x,;)  Ri  z(x)  +  {Vz(x)}T(xl  -  x) 

+  ^(x;  —  x)1  {'Hz(x)}(xi  —  x)  4 -  (2) 

=  do  +  di  (x;  -  x) 

+/3jvech{(xi  -  x)(xj  -  x)T}  -| - ,  (3) 

where  V  and  TL  are  the  gradient  and  Hessian  operators  respec¬ 
tively,  and  vech(-)  is  the  half-vectorization  operator  [2],  which 
lexicographically  orders  the  “lower-triangular”  portion  of  a  ma¬ 
trix  into  a  column  vector.  Indeed  the  local  approximation  can  be 
also  built  upon  bases  other  than  polynomials  [8] 

The  above  suggests  that  if  we  now  think  of  the  Taylor  se¬ 
ries  as  a  local  representation  of  the  regression  function,  estimat¬ 
ing  the  parameter  /30  can  yield  the  desired  (local)  estimate  of  the 
regression  function  based  on  the  data.  Indeed,  the  parameters 


{Pn}n= i  provide  localized  information  on  the  n-th  deriva¬ 

tives  of  the  regression  function.  Naturally,  since  this  approach  is 
based  on  local  approximations,  classical  regression  methods  es¬ 
timate  the  coefficients  {Pn}^=0  from  the  data  while  giving  the 
nearby  samples  higher  weights  than  samples  farther  away  (“ge¬ 
ometric”  weighting).  However,  it  is  also  appropriate  to  weight 
samples  based  on  their  relative  location  with  respect  to  a  local 
edge  (“radiometric”  weighting),  preforming  the  regression  along 
and  not  across  the  edges,  which  is  the  basis  of  modern  adaptive 
methods  .  A  general  formulation  we  propose,  capturing  this  idea 
is  to  solve  the  following  optimization  problem: 

p 

min  T,  m  ~  Po~  Pi  (xj  -  x)- 
/3^vech  |(Xj  -  x)  (xj  —  x)Tj - |  K(xj  -  x,  yt  -  y)  (4) 


where  K(-)  is  the  kernel  function  which  penalizes  both  geomet¬ 
ric  and  radiometric  distances  and  will  be  described  in  detail  in 
Section  3,  and  m  is  the  penalizing  parameter.  To  the  best  of  our 
knowledge,  all  kernel  regression  based  methods  in  the  literature 
choose  the  penalizing  parameter  as  m  =  2,  and  therefore  pose 
(4)  as  a  weighted  least-squares  problem.  In  Section  4,  we  show 
that  robustness  with  respect  to  the  outliers  can  be  significantly 
improved  by  exploiting  other  values  for  the  penalizing  parameter 
such  as  m  =  1,  which  in  effect  incorporates  a  robust  l\  norm  es¬ 
timator  [1]  in  the  kernel  regression  framework.  Furthermore,  we 
propose  novel  ways  to  adopt  the  kernel. 

Using  the  matrix  form,  the  optimization  problem  (4)  can  be 
posed  as  weighted  lm  norm: 

b  =  argmin||y-Xxb||™x,  (5) 


where 
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with  “diag”  defining  the  diagonal  elements  of  a  diagonal  matrix. 
We  use  steepest  descent  to  find  the  solution  to  this  minimization 
problem: 


,-("  +  !) 


b  =  b 


(») 


ffXj  Wxsign^y  -  Xxb<n))©|y  -  Xxb 


rwr 


(9) 


where  a  is  a  scalar  defining  the  step  size  in  the  direction  of  the 
gradient,  and  0  is  the  element  by  element  multiplication  operator. 

The  order  (N)  of  regression  affects  the  complexity  of  the 
local  approximation  of  the  signal.  In  the  non-parametric  statis¬ 
tics  literature,  locally  constant,  linear  and  quadratic  approxima¬ 
tions  (corresponding  to  N  =  0, 1,  2  respectively)  have  been  most 
widely  considered  [2],  In  particular,  choosing  local  constant  es¬ 
timation  with  m  =  2,  a  locally  linear  adaptive  filter  is  obtained. 


which  is  known  as  the  Nadaraya-Watson  Estimator  (NWE)  [3]. 
In  general,  lower  order  approximates,  such  as  NWE,  result  in 
smoother  images  (large  bias  and  small  variance)  as  there  are  fewer 
degrees  of  freedom.  On  the  other  hand  over-fitting  happens  in  re¬ 
gressions  using  higher  orders  of  approximation,  resulting  in  small 
bias  and  larger  estimation  variance.  Note  that,  in  the  experiments 
of  Section  4  we  used  the  second  order  (N  =  2)  approximation. 

3.  KERNEL  FUNCTION  SELECTION 

The  choice  of  kernel  function  greatly  affects  the  quality  of  re¬ 
construction.  In  this  section,  first  we  briefly  review  the  classic 
“non-adaptive”  kernel  function,  and  then  generalize  it  to  derive 
two  adaptive  kernel  functions  with  superior  performance. 

3.1.  Classic  Kernel  Function 

In  classic  kernel  regression,  samples  are  weighted  based  only  on 
their  spatial  distances  to  a  sample  of  interest,  which  simplifies  the 
kernel  K(-)  in  (4)  to 

K(xj  -  x,  i a  -  y )  =  Ah,  (x,;  -  x),  (10) 

where  A’h,  (•)  is  defined  as 

= skA'  K*)  •  <"• 

which  penalizes  distance  away  from  the  local  position  where  the 
approximation  is  centered.  The  2x2  “smoothing”  matrix  H, 
controls  the  strength  of  this  penalty.  The  standard  choice  of  the 
smoothing  matrix  is  Hi  =  hut I2,  where  pi  is  a  scalar  that  cap¬ 
tures  the  local  density  of  data  samples  and  h  is  the  global  smooth¬ 
ing  parameter,  extending  the  kernel  to  contain  “enough”  samples. 
As  described  in  [3],  in  case  of  irregularly  sampled  data,  it  is  rea¬ 
sonable  to  use  smaller  kernels  in  the  areas  with  more  available 
samples,  whereas  larger  kernels  are  more  suitable  for  the  more 
sparsely  sampled  areas  of  the  image.  The  choice  of  the  particular 
form  of  the  function  AT(-)  is  open,  and  may  be  selected  as  any 
symmetric  function,  which  attains  its  maximum  at  zero  such  as 
Gaussian. 

Since  the  shape  of  the  classic  kernels  is  independent  of  the 
radiometric  (gray  level)  information,  as  described  in  [3],  classic 
kernel  based  regression  methods  suffer  from  an  inherent  limita¬ 
tion  due  to  the  local  linear  action  on  the  data.  In  what  follows, 
we  discuss  extensions  of  the  kernel  regression  method  that  en¬ 
able  this  structure  to  have  nonlinear,  more  effective,  action  on  the 
data.  The  proposed  adaptive  kernel  functions  rely  on  not  only 
the  sample  location  and  density,  but  also  the  radiometric  proper¬ 
ties  of  these  samples.  Therefore,  the  effective  size  and  shape  of 
the  regression  kernel  are  adapted  locally  to  image  features  such 
as  edges.  This  property  is  illustrated  in  Fig.  1 .  where  the  clas¬ 
sical  and  adaptive  kernel  shapes  in  the  presence  of  an  edge  are 
compared. 

3.2.  Bilateral  Kernel  Function 

A  simple  and  intuitive  choice  of  the  adaptive  kernel  K(-)  is  to 
use  separate  terms  for  penalizing  the  spatial  and  radiometric  dis¬ 
tances.  Indeed  this  is  precisely  the  thinking  behind  the  bilateral 
filter,  introduced  in  [5,  6].  The  bilateral  kernel  choice  is  then 

K (x,  -x,yi-y)  =  A'Hi(x2  -  x)Khr(yi  -  y),  (12) 


Orientation  vector 


Fig.  1.  Kernel  spread  in  a  uniformly  sampled  data  set.  (a)  Kernels 
in  the  classic  method  depend  only  on  the  sample  density,  (b) 
Adaptive  kernels  elongate  with  respect  to  the  edge. 


would  otherwise  have  resulted,  is  tempered  around  edges  with 
appropriate  choice  of  C;’s.  With  such  steering  matrices,  for  ex¬ 
ample,  if  we  choose  a  Gaussian  kernel,  the  steering  kernel  is 
mathematically  represented  as 


Kh=  (xj  -  x) 


v/det(Ci) 

2nh2 


exp 


(x,  -  x)TCi(xj  -  x) 
2h2 


.  (15) 


The  local  edge  structure  is  related  to  the  gradient  covariance  (or 
equivalently,  the  locally  dominant  orientation).  In  [3]  we  have 
shown  that  a  convenient  form  of  representing  the  covariance  ma¬ 
trix,  is  to  decompose  it  into  three  components  as  follows: 


where  hr  is  the  radiometric  smoothing  scalar  that  controls  the 
rate  of  decay,  and  (•)  and  Kf,r{-)  are  the  spatial  and  radio- 
metric  kernel  functions,  respectively.  In  general,  the  application 
of  bilateral  kernel  is  limited  to  denoising  problem,  since  the  pixel 
value  (y)  at  an  arbitral  position  (x)  might  not  be  available  from 
data.  This  limitation,  however,  can  be  overcome  by  using  an  ini¬ 
tial  estimate  of  y  by  an  appropriate  interpolation  technique  [3]. 
Also,  breaking  K(-)  into  spatial  and  radiometric  terms  as  utilized 
in  the  bilateral  case  weakens  the  estimator  performance  since  it 
limits  the  degrees  of  freedom  and  ignores  correlations  between 
positions  of  the  pixels  and  their  values.  The  following  section 
provides  a  solution  to  overcome  this  drawback. 

3.3.  Steering  Kernel  Function 

Based  upon  the  earlier  non-parametric  framework,  the  filtering 
procedure  we  propose  next  takes  the  above  ideas  one  step  further. 
In  particular,  we  observe  that  the  effect  of  computing  Khr{yi  —  y) 
in  (12)  is  to  implicitly  measure  a  function  of  the  local  gradient  es¬ 
timated  between  neighboring  values,  and  to  use  this  estimate  to 
weight  the  respective  measurements.  As  an  example,  if  a  pixel  is 
located  near  an  edge,  then  pixels  on  the  same  side  of  the  edge  will 
have  much  stronger  influence  in  the  filtering.  With  this  intuition 
in  mind,  we  propose  a  two-step  approach  where  first  an  initial 
estimate  of  the  image  gradients  is  made  using  some  kind  of  gra¬ 
dient  estimator  (say  the  second  order  classic  kernel  regression 
method).  Next,  this  estimate  is  used  to  measure  the  dominant  ori¬ 
entation  of  the  local  gradients  in  the  image.  In  a  second  filtering 
stage,  this  orientation  information  is  used  to  adaptively  “steer” 
the  local  kernel,  resulting  in  elongated,  elliptical  contours  spread 
along  the  directions  of  the  local  edge  structure.  With  these  locally 
adapted  kernels,  the  denoising  is  effected  most  strongly  along  the 
edges,  rather  than  across  them,  resulting  in  strong  preservation  of 
details  in  the  final  output.  To  be  more  specific,  the  steering  kernel 
takes  the  form 

K(xj  —x.,yi  —  y)  =  (x,  —  x),  (13) 

where  Hips  are  the  data-dependent  full  matrices  which  we  call 
steering  matrices.  They  are  defined  as 

(14) 

where  CPs  are  (symmetric)  covariance  matrices  based  on  the  lo¬ 
cal  gray-values.  A  good  choice  for  C*  will  effectively  spread  the 
kernel  function  along  the  local  edges  as  shown  in  Fig.  1 .  It  is 
worth  noting  that  even  if  we  choose  a  large  h  in  order  to  have 
a  strong  denoising  effect,  the  undesirable  blurring  effect  which 


Ci=7<Ue4.Mj£,  (16) 
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where  Ugi  is  the  rotation  matrix  and  A;  is  the  elongation  matrix. 
Now  the  covariance  matrix  is  given  by  the  three  parameters  7 j, 
9i  and  ay,  which  are  the  scaling,  rotation,  and  elongation  param¬ 
eters,  respectively  and  the  effect  of  which  are  as  follows.  First, 
the  initial  circular  kernel  is  elongated  by  the  elongation  matrix 
A i  with  semi-minor  and  major  axes  given  by  oy  and  a~x ,  respec¬ 
tively.  Second,  the  elongated  kernel  is  rotated  by  the  matrix  U^. 
Finally,  the  kernel  is  scaled  by  the  scaling  parameter  7;.  We  re¬ 
fer  the  reader  to  [3]  for  the  details  of  estimating  these  parameters 
in  an  iterative  fashion.  We  note  that  the  presented  formulation  is 
close  to  the  apparently  independently  derived  normalized  convo¬ 
lution  formulation  of  [4], 

4.  EXPERIMENTS 

In  this  section  we  compare  the  performance  of  the  proposed  algo¬ 
rithm  to  other  methods.  We  show  that  in  presence  of  white  Gaus¬ 
sian  noise  the  proposed  robust  kernel  regression  method  works 
as  well  if  not  better  than  the  state  of  the  art  recent  wavelet  based 
denoising  method  of  [9],  and  other  popular  methods.  We  also 
note  that  the  wavelet  method  in  general  is  computationally  more 
efficient  than  the  steering  kernel  method.  However,  in  presence 
of  other  noise  models  (such  as  salt  and  pepper  noise)  while  the 
performance  of  non-robust  methods  dramatically  degrades,  the 
proposed  (1  based  robust  method  effectively  removes  the  noise. 
The  criterion  for  parameter  selections  in  all  the  examples  was  to 
choose  parameters  which  gave  the  best  RMSE  values. 

In  the  first  experiment,  we  added  white  Gaussian  noise  with 
standard  deviation  (STD)  of  25  to  the  original  image  of  Fig.  2(a) 
resulting  in  the  noisy  image  of  Fig.  2(b).  Denoised  images  us¬ 
ing  the  wavelet1  method  of  [9];  classic  kernel  regression  method 
{m  =  2,  h  =  1.33),  steering  kernel  regression  method  ( m  =  2, 
h ,  =  1.33,  7  iterations  initialized  with  Z2  classic),  steering  kernel 
regression  method  (m  =  1,  h  =  3,  2  iterations  initialized  with 
/]  classic)  and  corresponding  Root  Mean  Square  Error  (RMSE) 
values  are  shown  in  Fig.  2(c)-(f),  respectively. 

In  the  second  experiment  we  added  20%  salt  and  pepper  noise 
to  the  original  image  of  Fig.  2(a)  resulting  in  the  noisy  im¬ 
age  of  Fig.  3(a).  Denoised  images  using  a  3  x  3  median  filter, 
wavelet  method  of  [9],  classic  kernel  regression  method  (m  =  2, 

This  result  is  produced  by  the  software,  available  on 
http://decsai.ugr.es/~javier/denoise/index.html. 


(d)  l2  Classic  (e)  l2  Steering  (f)  l\  Steering 

Fig.  2.  Gaussian  noise  removal  experiment.  Corresponding 
RMSE  values  for  (b)-(f)  are  25.0,  9.71,  11.36,  10.11,  and  10.71, 
respectively. 


Fig.  3.  Salt  &  pepper  noise  removal  experiment.  Corresponding 
RMSE  values  for  Figures(a)-(f)  are  63.84,  11.05,  21.54,  21.81, 
21.06,  and  7.14,  respectively. 

h  =  2.46),  steering  kernel  regression  method  (m  =  2,  h  =  2.25, 
20  iterations  initialized  with  l2  classic),  steering  kernel  regression 
method  (m  =  l,h  =  2.25,  zero  iteration  initialized  with  l\  clas¬ 
sic)  and  corresponding  RMSE  values  are  shown  in  Fig.  3(b)-(f), 
respectively. 

In  our  final  experiment,  we  added  white  Gaussian  noise  with 
STD  of  10  along  with  5%  salt  and  pepper  noise  to  the  original 
image  of  Fig.  2(a).  Then,  we  randomly  eliminated  85%  of  these 
noisy  pixels,  creating  the  sparse  image  of  Fig.  4(a).  Interpo¬ 
lated  and  denoised  images  using  the  Delaunay-spline  smoother 
(refer  to  [3]  for  details),  and  the  iterative  steering  kernel  regres¬ 
sion  method  (m  =  1,  h  =  3,  0  iterations)  and  corresponding 
RMSE  values  are  shown  in  Fig.  4(b)-(c),  respectively. 

5.  CONCLUSIONS 

In  this  paper  we  promoted,  extended,  and  demonstrated  kernel 
regression  as  a  general  framework  for  studying  several  efficient 
denoising  and  interpolation  algorithms.  To  overcome  the  inherent 
limitations  dictated  by  the  linear  filtering  properties  of  the  clas¬ 
sic  kernel  regression  methods,  we  introduced  the  non-linear  data- 
adapted  class  of  kernel  regressors  with  superior  performance.  Fur¬ 


(a)  Subsampled  (b)  Del.  Spline  (c)  Zj  Steering 


Fig.  4.  Sparse-noisy  image  interpolation  experiment,  (a)  is  the 
input  image  with  85%  of  pixels  removed,  and  further  corrupted 
by  adding  Gaussian  and  salt  and  pepper  noise.  Reconstructed  im¬ 
ages  using  the  Delaunay-spline  smoother  (RMSE=22.5),  and  the 
li  steering  kernel  regression  (RMSE=17.5)  methods,  are  shown 
in  (b)-(c),  respectively. 

thermore,  we  achieved  robustness  with  respect  to  outliers  in  data 
and  noise  model  by  incorporating  the  Zj  norm  penalty  in  the  ker¬ 
nel  regression  framework.  Image  deblurring  is  also  an  important 
issue  in  image  reconstruction,  and  it  is  a  part  of  our  ongoing  work 
within  this  framework. 
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ABSTRACT 

We  exploit  recent  advances  in  the  physical  design  of  fast  op¬ 
tical  systems  which  enable  active  imaging  with  “ballistic” 
light.  In  this  modality,  fast  bursts  of  optical  energy  are  prop¬ 
agated  into  a  medium,  and  the  ballistic  component  of  light 
(which  travels  with  minimal  diffusive  distortion)  is  detected 
after  transmission  through  the  target  and  the  medium.  To  im¬ 
prove  the  detection  rate  of  the  common  single  pixel  optimal 
detectors,  we  exploit  sampling  at  a  diversity  of  locations  in 
space,  and  develop  a  multi-scale  algorithm  based  upon  the 
Generalized  Likelihood  Ratio  Test  (GLRT)  framework,  which 
takes  advantage  of  the  spatial  correlation  of  nearby  samples. 
Experimental  results  show  that  objects  of  different  size  and 
shape  that  are  completely  unrecognizable  using  the  common 
single  pixel  detection  techniques,  are  detectable  with  very 
high  accuracy  using  the  said  multi-scale  GLRT  technique. 

Index  Terms —  Ballistic  Photons,  Poisson  Statistics,  Adap¬ 
tive  Reconstruction,  GLRT,  Coherent  Imaging,  Turbid  Media. 

1.  INTRODUCTION 

High  resolution  imaging  and  detection  of  objects  hidden  in 
a  turbid  (scattering)  medium  have  long  been  challenging  and 
important  problems  with  many  industrial,  military,  and  med¬ 
ical  applications.  While  turbid  media  such  as  fog,  smoke, 
haze,  or  body  tissue  are  virtually  transparent  to  radar  range 
electromagnetic  waves,  the  resolution  of  radar-based  imaging 
systems  is  often  insufficient  for  many  practical  applications. 
On  the  other  hand,  while  the  resolution  of  imaging  systems 
using  ultra  short  wavelengths  (e.g.  X-rays)  is  very  desirable, 
there  exist  potential  health  hazards  for  imaging  subjects  and 
technicians  alike. 

As  an  alternative,  imaging  systems  working  in  the  optical/ 
infra-red  spectrum  range  (laser  scanners)  are  potentially  able 
to  produce  high  resolution  images  without  the  likely  health 
hazards.  Unfortunately,  even  a  very  thin  and  powerful  colli¬ 
mated  laser  beam  quickly  diffuses  as  it  travels  in  turbid  me- 
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dia,  similar  to  a  car’s  headlights  in  fog.  Therefore,  a  naive 
approach  to  optical  imaging  of  objects  hidden  inside  a  turbid 
medium  results  in  very  blurry  images,  where  targets  are  often 
indistinguishable  from  each  other  or  the  background. 

Fortunately,  the  advent  of  the  new  tunable  solid  state  lasers 
and  ultra-fast  optical  detectors  has  enabled  us  to  acquire  high- 
quality  images  through  turbid  media  where  the  resolution  is 
only  limited  by  diffraction.  While  many  efficient  imaging 
systems  for  capturing  high-resolution  images  through  turbid 
media  have  been  proposed  through  the  years,  in  this  paper 
we  mainly  focus  on  ultra-fast  time-gated  or  coherent  imaging 
systems  first  introduced  in  [  1  ] . 

Ultra-fast  time-gated  imaging  is  based  on  scanning  the 
region  of  interest  (ROI)  point-by-point  by  sending  very  fast 
bursts  of  optical  energy  (laser  pulses)  and  detecting  the  un¬ 
scattered  (coherent)  photons  that  have  passed  through  the 
medium  or  reflected  from  the  object.  Although  most  of  the 
photons  in  a  laser  pulse  are  either  randomly  scattered  (los¬ 
ing  their  coherence)  or  absorbed  as  they  travel  through  tur¬ 
bid  media,  across  short  distances,  a  few  photons  keep  their 
coherence  and  pass  through  in  straight  lines  without  being 
scattered.  These  coherent  photons  are  commonly  referred  to 
as  the  ballistic  photons.  Aside  from  the  diffusive  and  ballistic 
photons,  the  photons  that  are  slightly  scattered  retaining  some 
degree  of  coherence  are  referred  to  as  snake  photons.  Since 
the  diffusive  and  ballistic  photons  have  different  path  lengths, 
a  femto-second  laser  pulse  generator  and  an  ultra  fast  time 
gate  can  be  paired  to  separate  the  relatively  slow  (delayed) 
diffusive  photons  from  the  ballistic  ones. 

In  what  follows  in  this  paper,  we  focus  on  studying  and 
improving  the  performance  of  ballistic  imaging  systems.  In 
Section  2,  we  describe  a  statistical  model  for  the  signal  and 
noise  in  a  typical  ballistic  imaging  scenario.  In  Section  3,  we 
study  optimal  single  pixel  detection  systems  and  show  that 
better  detection  rates  are  achievable  using  a  multi-pixel  detec¬ 
tion  technique  which  is  based  on  the  GLRT  principle.  Section 
4  concludes  this  paper. 

2.  STATISTICAL  MODEL  FOR  BALLISTIC 
IMAGING  SYSTEMS 

To  have  a  better  understanding  of  the  practical  issues  involved 
in  photon  limited  imaging  via  ballistic  systems,  let  us  con- 


sider  the  imaging  system  described  in  [2],  where  the  pumped 
Ti:Sapphire  laser  radiates  800nm  pulses  at  a  repetition  rate 
of  1  kHz  and  an  average  power  of  60m W.  It  is  easy  to  show 
that  the  number  of  photons  in  each  packet  of  energy  (pulse)  is 
computed  as 


Pulse  Energy 
Photon  Energy 


60x  10~3  x  Is 
_ 1000 _ 

2.4830  x  10~19 


=  2.4164xl014.  (1) 


Due  to  the  statistical  nature  of  pulse  propagation,  as  a 
laser  beam  travels  through  a  diffusive  medium,  it  is  possi¬ 
ble  that  some  of  the  photons  emerge  without  being  scattered. 
By  selecting  these  unscattered  “ballistic”  photons,  and  reject¬ 
ing  the  scattered  (diffused)  ones,  it  is  possible  to  obtain  non- 
blurred  images  which  are  the  sharp  shadows  of  targets  buried 
in  the  diffusive  medium. 

As  expected,  in  relatively  long  distances,  the  number  of 
detected  ballistic  photons  is  extremely  small.  Indeed,  Beer’s 
Law  dictates  an  exponential  relationship  between  the  intensity 
of  the  transmitted  light,  and  that  of  the  ballistic  component  as 

h  =  /0exp(-^).  (2) 

In  this  expression,  Iq  is  the  number  of  the  generated  photons 
in  one  laser  pulse  before  entering  the  turbid  medium,  Ib  is 
the  number  of  the  ballistic  photons  which  survive  traveling 
through  the  medium,  d  is  the  distance  traveled  through  the 
medium,  L  =  —  is  the  mean  free  path  (MFP)  length  (average 
distance  photons  travel  before  being  scattered),  and  jit  is  the 
medium  extinction  factor.  From  (1)  and  (2),  it  is  clear  that 
for  typical  laser  powers,  it  is  fairly  unlikely  that  any  ballistic 
photon  survive  imaging  scenarios  where  the  ratio  of  d/L  is 
larger  than  ~  30  MFP. 

The  exponential  drop  in  the  number  of  received  photons 
is  the  main  prohibitive  factor  for  using  such  high-resolution 
optical  imaging  systems  across  long  distances.  In  such  imag¬ 
ing  scenarios,  we  are  forced  to  rely  on  the  less  informative 
snake  and  diffusive  photons.  In  [3],  an  accurate  yet  compu¬ 
tationally  manageable  mathematical  model  for  diffusive  light 
propagation  in  turbid  media  is  presented.  An  example  of  such 
imaging  modality  and  experimental  analysis  is  presented  in 
[4]  and  some  excellent  literature  surveys  on  the  subject  of  dif¬ 
fusive  imaging  systems  are  presented  in  [5].  However,  imag¬ 
ing  systems  that  are  able  to  time-resolve  both  ballistic  and 
diffusive  photons  are  rather  expensive  and  are  not  discussed 
in  this  paper.  Here,  we  focus  on  imaging  systems  that  detect 
ballistic  photons  only.  We  exploit  these  statistical  studies  to 
improve  the  performance  of  ballistic  imaging  systems  even  in 
long  distances  where  the  signal  power  is  weak. 

It  is  important  to  note  that  due  to  the  stochastic  nature  of 
photon  propagation,  Ib,  calculated  in  (2),  is  merely  the  ex¬ 
pected  value  of  a  Poisson  random  variable  that  estimates  the 
number  of  surviving  ballistic  photons.  Moreover,  we  assume 
that  the  received  signal  at  the  detector  is  contaminated  with 
some  amount  of  independent  Poisson  noise  due  to  shot  noise 


and  other  degrading  effects.  Therefore,  since  the  received  sig¬ 
nal  at  the  detector  is  the  unweighted  summation  of  two  Pois¬ 
son  random  variables,  it  can  be  modeled  as  a  Poisson  random 
process  with  the  following  expected  value 

/  =  /o  e  +  Xe  =  Xs  +  Xe, 


where  Xe  and  Xs  are  the  expected  values  of  the  noise  and 
signal,  respectively. 

Considering  such  imaging  model,  the  probability  density 
function  of  the  received  signal  is  given  by 
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f(y\X^  +  X.)=]J 


k= 1 


e~^+x^(Xek+XSk) 

Vk]- 


Vk 


(3) 


where  yk  is  the  k-th  measurement,  y=  [j/i,  y2,  yk,  — ,  Vn}T, 
X c  [-^ei  7  Xe2 , ...,  Xek , ...,  XeN]  ,  and 

Xs  =  [XSl,XS2,  ...,XSk, ...,XSn]t.  Note  that  the  laser  emits 
thousands  of  pulses  per  second  and  in  practical  implementa¬ 
tion  each  spatial  position  is  measured  N  times  to  improve  the 
quality  of  estimation,  and  therefore  the  model  in  (3)  is  pre¬ 
sented  in  the  vector  form.  Since  the  average  power  of  laser  or 
the  detector  (and  medium)  characteristics  are  assumed  not  to 
be  changing  abruptly,  to  simplify  notations,  we  assume  that 
Xei  —  Xe2  —  ...  —  XeN  —  Xe ,  and  XSl  —  XS2  —  ...  — 
XSN  =  Xs  (extension  to  the  more  general  time-varying  sig¬ 
nal  and  noise  case  is  straight  forward). 

3.  OPTIMAL  DETECTION  OF  OPAQUE  OBJECTS 
IN  TURBID  MEDIA 


In  this  section,  assuming  that  the  medium,  laser,  target,  and 
turbid  medium  are  accurately  calibrated,  we  present  the  sta¬ 
tistical  optimal  detectors  of  opaque  objects  hidden  in  a  turbid 
medium. 

3.1.  Single  Pixel  Optimal  Detection 

In  this  subsection,  we  study  the  Neyman-Pearson  (N.P.)  type 
statistical  test  [6]  for  detecting  opaque  objects  hidden  in  a  tur¬ 
bid  medium.  In  this  test,  we  basically  compare  the  likelihood 
of  the  following  two  scenarios: 

•  H0 :  An  opaque  object  is  hidden  in  the  scattering  medium, 
blocking  the  laser  pulse  (  i.e.  measurements  contain 
only  noise). 

•  Hi :  No  opaque  object  exists  in  the  propagation  line  of 
the  laser  pulse  (i.e.  measurements  contain  noise  plus 
attenuated  laser  pulse). 


The  probability  density  function  of  these  two  scenarios  when 
such  tests  are  repeated  N  times  are  given  by 

Ho  :  f(y\Xe)  =  l[-  1  e) 
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and  therefore  the  N.P.  test  is  derived  by  comparing  the  log 
likelihood  ratio  to  a  threshold  as: 
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Noting  that  ^  is  yet  another  poisson  process,  the  proba- 

k=l 

bilities  of  false  alarm  (Pfa)  and  detection  (Pd)  are  computed 
as 

N  y  p-NXe(  ATV  \k 
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Fig.  1.  An  illustrative  example  showing  the  outline  of  the  pro¬ 
posed  multi-scale  GLRT  algorithm,  a:  Scale  1  measurement. 
b,c:  Scales  2,3  super-pixels,  respectively,  d-f:  Confidence 
values  for  scales  1-3,  respectively.  The  check  marked  second 
scale  gives  the  highest  confidence  value  for  the  central  pixel. 


3.2.  Multi-pixel  GLRT  Detection 

As  explained  in  Section  1,  in  ballistic  imaging  the  field  of 
view  is  scanned  at  multiple  points  to  create  a  2-D  image  of 
the  objects  in  the  ROI.  In  this  section,  we  propose  an  effective 
algorithm  that  exploits  the  spatial  correlation  of  the  nearby 
samples  in  a  multi-pixel  imaging  scenario  to  improve  on  the 
performance  of  the  single  pixel  optimal  detectors  developed 
in  the  previous  section. 

The  proposed  multi-pixel  detection  technique  generalizes 
the  single  pixel  detection  techniques  and  preforms  optimal 
tests  on  “super-pixels”,  which  are  the  collective  intensities  of 
a  set  of  neighboring  pixels  in  size  and  shape  of  the  hidden 
objects.  However,  since  in  general  the  size  and  shape  of  the 
hidden  objects  is  not  known  a  priori,  we  develop  a  GLRT 
based  algorithm  that  simultaneously  tests  the  existence,  and 
estimates  the  shape  and  size  of  the  objects  hidden  in  turbid 
media. 

The  outline  of  the  proposed  GLRT  algorithm  is  illustrated 
by  an  example  in  Fig.l.  First,  for  a  given  (fixed)  false  alarm 
rate  the  optimal  detectors  developed  in  the  previous  section 
are  exploited  to  test  the  existence  or  absence  of  objects  at  each 
individual  pixel.  As  an  illustrative  example,  this  test  is  applied 
to  the  central  pixel  (shaded)  of  Fig.  1(a),  where  the  measured 
pixel  value  (0.4)  is  compared  to  the  N.R  test  threshold  (0.5). 
Of  course,  the  greater  the  distance  of  the  measurement  from 
the  threshold,  the  more  confident  we  are  in  the  accuracy  of 
the  test  result.  Next,  we  integrate  the  gray-level  values  of  all 
immediate  neighboring  pixels,  and  in  effect  consider  them  as 
one  “super-pixel”,  as  illustrated  in  Fig.  1(b).  Since  the  false 
alarm  rate  is  fixed  for  all  scales,  the  decision  threshold  is 
different  than  the  threshold  calculated  in  the  previous  step, 
which  is  recalculated  based  on  the  grayvalue  of  the  super¬ 
pixel.  In  the  next  steps,  we  repeat  this  process  by  fixing  the 


false  alarm  rate  and  considering  larger  neighborhoods.  The 
generalized  N.R  test  for  these  steps  is  formulated  as  follows 


scale 

Vm,l 


H,  log(7«cate)  +  NgcaleXs 

Ho  iog(7f7) 
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where  y^°ile  is  the  summation  of  the  pixel  values  in  the 
N scale  =  N( 2x  scale—  l)2  pixels  neighborhood  around  the 
pixel  [to,  l]  (note  that  other  neighborhood  expansion  strate¬ 
gies  with  different  shape  and  size  can  be  also  considered  in 
this  algorithm).  Our  confidence  in  the  decision  made  on  each 
scale  is  defined  as  the  distance  between  the  summation  of 
measurements  in  the  super-pixel  and  that  of  the  threshold  set 
by  the  GLRT: 


Confidence^//6 


I V 


scale 

m,l 


log(7scaJ  +  NscaleXs 


log( 


Xe+Xs 

Xe 


) 


(9) 


Note  that  the  optimal  scale  is  not  unique  for  all  pixels,  as  finer 
scales  are  more  suitable  for  pixels  located  on  the  texture  or 
edge  areas,  and  coarser  scales  are  more  suitable  for  the  pixels 
located  in  flat  areas.  Therefore,  we  decide  on  the  presence  or 
absence  of  the  object  at  a  particular  pixel  based  on  the  test  re¬ 
sult  of  the  scale  that  shows  the  highest  confidence  value.  The 
memory  requirement  of  this  technique  is  independent  of  the 
maximum  scale  number,  since  we  only  need  to  keep  the  orig¬ 
inal  image,  the  last  estimated  image  and  the  corresponding 
confidence  values. 

To  have  a  better  understanding  of  the  proposed  multi-scale 
GLRT  technique  and  its  performance,  we  set  up  an  illustra¬ 
tive  controlled  imaging  scenario.  Fig. 2(a)  shows  ideal  (noise¬ 
less  and  deterministic)  image  of  objects  of  different  size  and 
shape.  To  depict  an  experiment  at  the  limit  distance  where  the 
signal  of  interest  is  very  weak,  we  consider  an  imaging  sce¬ 
nario  where  the  average  number  of  received  ballistic  photons 


for  each  pixel  is  one  photon.  Fig. 2(b)  shows  such  Poissonian 
random  signals  (yet  free  of  noise  effects).  Detection  of  such 
signals  becomes  more  difficult  when  we  consider  the  system 
noise  as  illustrated  in  Fig.2(c),  where  the  Poisson  noise  vari¬ 
ance  is  40.  Fig.2(d)  is  an  image  reconstructed  by  implement¬ 
ing  the  point-by-point  single  pixel  detection  techniques,  con¬ 
sidering  a  false  alarm  rate  of  0.00125,  where  none  of  the  ob¬ 
jects  are  correctly  identified.  On  the  other  hand.  Fig. 2(e)  is  the 
result  of  exploiting  the  multi-scale  GLRT  technique,  showing 
a  considerably  more  accurate  detection  of  such  objects. 

Fig. 2(f)  illustrates  the  scale  from  which  each  pixel  in  the 
final  image  of  Fig. 2(e)  is  selected.  Note  that  as  expected,  the 
pixels  in  the  flat  area  are  selected  from  the  coarser  scales, 
while  the  pixels  on  the  edge  areas  are  selected  from  the  finer 
scales.  Fig. 2(g)  shows  the  confidence  in  the  detection  re¬ 
sult  (9)  with  respect  to  the  corresponding  pixels.  This  figure 
shows  higher  confidence  levels  in  the  flat,  and  less  confidence 
in  the  edge  areas.  Also,  in  Fig. 2(g)  we  see  that  the  area  with 
the  lowest  confidence  is  the  place  where  most  misclassifica- 
tions  happen.  This  is  good  news,  since  to  increase  the  detec¬ 
tion  rate,  we  may  opt  to  do  a  second  (and  very  quick)  round  of 
scans,  sampling  only  on  these  very  low-confidence  regions.  In 
Fig. 2(h),  we  plot  the  misclassification  rates  at  each  scale  (blue 
line),  and  compare  it  to  the  overall  multi-scale  one  (red  line). 
These  experimental  plots  show  that  the  performance  of  the 
proposed  pixelwise  GLRT  technique  (depending  on  the  noise 
level)  is  either  very  close  or  even  better  than  the  best  fixed 
scale  technique.  In  Fig.2(i),  the  performance  of  single  pixel 
detection  technique  is  compared  with  the  multi-scale  ones  via 
their  corresponding  ROC  curves  (with  25  Monte  Carlo  exper¬ 
iments).  Once  again,  the  multi-scale  technique  shows  the  best 
or  close  to  the  best  performance. 

4.  CONCLUSION 

In  this  paper,  we  studied  a  technique  for  improving  the  qual¬ 
ity  of  the  ballistic  images  captured  through  turbid  media.The 
novelty  of  this  paper  is  in  combining  the  recent  advances  in 
optical  science  with  the  novel  image  processing  and  statistical 
signal  processing  techniques. 
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ABSTRACT 

The  tradeoffs  between  ballistic  imaging  (time-gated  imaging 
of  first-arrival,  unscattered  photons)  and  conventional  imag¬ 
ing  for  resolving  tumors  in  biological  scattering  media  are 
examined.  For  ballistic  imaging,  closed  form  expressions  are 
derived  to  characterize  the  resolvability  using  five  degrees  of 
freedom  (laser  intensity,  scattering  coefficient,  thickness  of 
medium,  false  alarm  rate,  and  number  of  observations).  For 
conventional  imaging,  a  numerical  approximation  is  used  to 
find  the  asymptotic  resolution  using  the  scattering  and  ab¬ 
sorption  coefficients  of  the  medium.  Using  the  characteri¬ 
zations  of  both  approaches,  a  decision-theoretic  approach  to 
determining  the  minimum  resolvable  object  size  is  developed, 
which  provides  clear  guidelines  as  to  when  time-gated  ballis¬ 
tic  imaging  methods  offer  advantages  over  conventional  imag¬ 
ing.  The  theoretical  predictions  are  validated  through  a  real¬ 
istic  simulation  of  tumors  in  breast  tissue. 

Index  Terms —  Transillumination  Imaging,  Decision  The¬ 
oretic,  Resolution,  Ballistic 

1.  INTRODUCTION 

Ballistic  photon  imaging  is  a  promising  methodology  for 
studying  highly  scattering  media  such  as  human  tissue.  Re¬ 
cent  advances  allow  for  the  time-gating  of  early  arriving  pho¬ 
tons  introduced  into  a  scattering  (turbid)  medium.  This  allows 
for  the  ability  to  separately  detect  unscattered  or  ballistic  pho¬ 
tons  that  exit  the  medium.  Due  to  the  lack  of  scattering,  these 
photons  will  retain  the  spatial  information  of  the  medium.  A 
tradeoff  occurs,  as  the  number  of  ballistic  photons  decays  ex¬ 
ponentially  fast  as  the  thickness/depth  of  the  turbid  medium 
increases.  This  results  in  an  observation  from  ballistic  pho¬ 
tons  that  offers  a  high  resolution  but  low  SNR. 

The  basic  imaging  model  considered  here  is  a  single  laser 
point  source  and  a  single  photon  detector  placed  on  either 
side  of  a  turbid  (scattering)  medium.  We  assume  that  the 
source  and  detector  can  be  positioned  at  arbitrary  points,  to 
allow  probing  through  any  desired  set  transects  through  the 
medium.  A  turbid  medium  can  be  considered  any  scatter¬ 
ing  material,  and  in  this  specific  analysis  we  assume  it  to  be 
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a  homogeneously  scattering  section  of  human  tissue.  This 
scattering  material  will  be  parameterized  by  the  number  of 
scatters  per  unit  length  (/tt™),  the  fraction  of  absorptions  per 
length  (/t,™),  and  the  length  of  the  medium  ( dm ).  At  certain 
points  between  the  source  and  detector  in  the  medium  there 
may  exist  an  occluder  of  interest.  The  occluders,  which  are 
assumed  to  be  opaque,  represent  tumors  embedded  in  tissue 
for  a  biomedical  imaging  application.  We  assume  that  the 
occluders  are  located  at  the  mid-point  of  the  medium  being 
probed,  as  this  is  the  point  in  the  medium  where  the  variance 
of  the  scattering  will  be  at  a  maximum  [1,2].  To  completely 
define  the  environment,  one  must  take  into  consideration  the 
scattering  properties  of  the  occluding  tumor  (//.[,,  //*),  and  the 
physical  properties  of  the  occluding  tumor  (depth  =  df ,  width 
=  w *). 


Turbid  Medium 


Source 

Point 

W\ 

Detector 

Point 

*~~v 

dm 

Figure  1.  Turbid  Medium  Model  Example 


To  create  an  observation  of  the  environment,  the  source 
and  detector  pair  will  make  a  raster  scan  of  the  field  of 
view  to  create  an  estimated  cross-section  of  the  environment. 
Our  goal  is  to  find,  given  the  parameters  of  the  medium 
and  the  occluder,  lower  bounds  on  the  width  and  depth  of 
the  occluder  ( dt ,  v/)  that  can  be  reliably  resolved  in  the 
medium.  Two  imaging  regimes  will  be  considered,  the  bal¬ 
listic  regime  where  only  early  arriving  photons  are  detected 
using  high-speed  gating  mechanisms,  and  the  conventional 
imaging  regime,  where  all  the  photons  arriving  at  the  detec¬ 
tor  are  utilized  (i.e.,  no  time-gating).  The  advantage  of  the 
ballistic  photons  is  that  they  are  not  scattered,  providing  very 
high  spatial  resolution.  The  limitation  of  ballistic  imaging 
is  that  very  few  photons  will  propagate  through  the  medium 
without  scattering,  resulting  in  a  very  poor  signal-to-noise  ra- 


tio  (SNR).  Conventional  imaging  uses  all  photons  (scattered 
and  unscattered)  and  thus  provides  a  complementary  trade-off 
(lower  resolution  but  higher  SNR).  A  decision-theoretic  tech¬ 
nique  is  then  derived  to  determine  which  technique  should  be 
used  for  a  given  parameterized  scattering  medium  in  order  to 
resolve  small  tumors. 

2.  BALLISTIC  ANALYSIS 

The  ballistic  imaging  regime  consists  of  probing  points  in  a 
Field-of- View  ( FOV )  and  acquiring  a  count  of  the  number  of 
ballistic  photons  that  arrive  at  each  point  in  the  FOV .  These 
photons  traverse  in  a  straight  line-of-sight  between  the  source 
point  and  the  detector,  and  therefore  retain  the  spatial  charac¬ 
teristics  of  the  occluding  objects  in  the  turbid  medium.  Due 
to  their  direct  line-of-sight  characteristics,  these  photons  will 
travel  the  shortest  length  and  will  arrive  at  the  detector  be¬ 
fore  any  scattering  photons.  In  order  to  collect  only  ballistic 
photons,  time-gating  is  performed  to  restrict  the  observation 
to  only  early  arriving  photons  (i.e.  photons  with  no  scattering 
events).  If  a  perfect  occluder  is  in  the  line-of-sight  between 
the  source  and  detector,  then  no  ballistic  photons  will  arrive 
at  the  detector.  Therefore,  the  detection  of  ballistic  photons 
indicates  a  “clear”  line-of-sight  and  the  absence  of  a  tumor 
along  the  given  transect.  However,  a  problem  occurs  -  as  the 
turbidity  of  the  medium  increases,  the  less  likely  a  photon  is 
to  arrive  at  the  sensor  having  no  scattering  events.  As  a  conse¬ 
quence,  the  total  number  of  ballistic  (non-scattering)  photons 
for  a  medium  might  be  very  low.  In  addition,  stray  “noise” 
photons  from  other  sources  corrupt  the  observed  signal,  re¬ 
sulting  in  an  observation  that  has  high  spatial  resolution,  but 
low  SNR. 

2.1.  Ballistic  Imaging  -  Single  Point 

The  detected  time-gated  photons  are  ballistic  photons  from 
the  source  or  stray  “noise”  photons  arriving  during  the  time- 
gate  interval  from  other  ambient  light  sources.  The  problem 
of  determining  whether  or  not  a  tumor  lies  along  the  line-of- 
sight  can  be  cast  as  a  statistical  hypothesis  test,  as  follows. 
Given  a  received  photon  count  X  at  the  sensor,  one  must 
choose  between  two  possible  situations.  The  first  situation, 
is  that  no  occluding  tumors  exist  in  the  path  between  the  laser 
and  the  sensor  (Ho)-  The  alternative  is  that  there  is  an  occlud¬ 
ing  tumor  along  the  line-of-sight  between  the  laser  and  the 
sensor  (Hi).  Tumors  can  be  considered  to  be  significantly 
more  turbid  scattering  medium  than  the  healthy  tissue  [3], 
and  therefore  the  detector  will  collect  an  attenuated  number 
of  ballistic  photons  (compared  to  Ho)  along  with  the  noise 
photons. 

We  define  the  number  of  noise  photons  that  will  arrive  at 
the  detector  as  V  (t\0)  (where  X  ~  V( A)  is  a  Poisson  dis¬ 
tributed  random  variable  with  mean  =  A,  and  where  t  is  the  du¬ 
ration  of  photon  acquisition  time-gating).  The  ballistic  pho¬ 


tons  traveling  through  only  healthy  tissue  will  be  V  {t\L(tissUe)) 
(where  A L(tissue)  =  Xl  (exp  (— prTOdm))),  while  the  ballistic 
photons  traveling  through  both  healthy  tissue  and  tumor  will 

be  V  (tXiJ(tiSSUejrturnor')'j  (where  A 

L(tissue-\-  tumor) 

X l  (exp  (—  —  nmdl  +  pt*d*))),  with  Xl  is  the  expected 

number  of  photons  sent  through  the  medium  by  the  laser  per 
unit  time,  pj im  =  pi™  +  pt™,  and  pi*  =  pt*  +  pt*  )• 

This  can  be  expressed  as  a  hypothesis  test  with  the  null  hy¬ 
pothesis  (only  tissue)  defined  as  Ho  :  X  ~  V  (t  (Ao  +  A L(tissue 
and  the  true  hypothesis  (tissue  and  tumor)  defined  as  H\  : 

X  ~  V(t(Xo  +  XL{tissue+tumor))) ■  As  the  mean  of  the 
Poisson  distribution  grows,  the  probability  distribution  tends 
to  a  Gaussian.  Averaging  repeated  trials  (i.e.  averaging  of 
multiple  laser  pulses)  therefore  results  in  a  Gaussian  distributed 
statistic.  Using  the  Anscombe  Transformation  [4],  we  obtain 
the  following  relationship  (where  X  ~  J\f  (pt,  a2)  is  a  gaus- 
sian  distributed  random  variable  with  mean  =  p;  and  variance 

=  cr2)  then  X  ~  V  (A)  =>•  2y/x+  §  ~  A f  (2y/\,  l) .  Defin¬ 
ing  a  new  variable  representing  the  Anscombe  Transformed 
statistic  X'  =  2^/x  +  |  the  hypothesis  test  becomes: 

Ho  :  X'  ~  X  (2  ft  (Ao  +  AL(tiMue)),l)  (D 
Hi  :  X'  ~  AT  ^2 \Jt  (Ao  +  Xi,(ti ssue+tumor))  i  (2) 

The  decision  test  is  now  defined  as  ^X'  V)  .  Using 
the  test,  a  user-specified  false  alarm  rate  (a)  determines  the 
value  of  the  threshold  (7*)  such  that  P  (X*  <  7*  |  Ho)  <  a. 

2.2.  Ballistic  Imaging  -  K  Points 

The  problem  now  is  modified  to  trying  to  image  a  fixed  square 
array  of  (\fK  x  \/~K)  points.  This  results  in  a  multiple  hy¬ 
pothesis  testing  problem  (K  tests),  and  for  large  K  it  is  diffi¬ 
cult  to  control  the  overall  probability  of  false-alarm  (i.e.,  false 
tumor  detection  at  one  or  more  point).  A  standard  technique 
is  to  increase  the  acquisition  time  (t)  for  each  point,  but  by 
the  setup  of  the  problem  increasing  t  will  increase  not  only 
the  number  of  signal  photons,  but  also  the  number  of  noise 
photons.  This  puts  a  lower  bound  on  the  SNR  of  the  obser¬ 
vation.  To  boost  the  SNR,  one  could  use  spatial  aggregation 
by  averaging  over  a  number  of  observation  points.  This  mod¬ 
ifies  the  problem  to  averaging  neighborhoods  of  points  in  an 
area  measuring  \fM  x  \J~M ,  M  <  K,  effectively  reducing 
the  spatial  resolution  of  the  detection  map  (image).  By  de¬ 
creasing  the  spatial  resolution,  this  also  decreases  the  variance 
at  each  point,  modifying  the  decision  test  to: 

Ho  :  X'  ~  J\f  ^2 \Jt  (Ao  +  AL(tissMe)),  — ^  (3) 

Hi  :  X'  ~  N  {^\jt  (Ao  +  A  L(tissue+tumor))  >  (4) 


This  test  is  under  the  assumption  that  the  averaging  win¬ 
dow  will  contain  either  no  occluder  points  or  all  occluder 
points.  In  reality,  the  averaging  filter  will  result  in  an  observed 

point:  X”  ~  Af(pE[X'\n0]  +  (l-p)E[X'\n1],^) 

(where  p  is  the  fraction  of  the  window  containing  non¬ 
occluders).  Our  goal  is  to  find  the  lower  bound  on  the  value  of 
M  that  will  guarantee  an  overall  false  alarm  rate  of  less  than 
a,  we  consider  the  ideal  case  (all  occluders  or  non-occluders) 
in  our  calculations  in  order  to  obtain  closed-form  solutions. 


2.2.1.  Bonferroni  Correction 


The  Bonferroni  Correction  approach  is  a  conservative  method 
of  controlling  the  false  alarm  rate  for  a  detection  problem  un¬ 
der  multiple  i.i.d.  tests  [5],  The  correction  adjusts  the  thresh¬ 
old  for  each  individual  test  in  order  to  satisfy  a  lower  (per  test) 
false  alarm  rate  value  (-^)  such  that  each  of  the  fixed  number 
K-points  in  the  array  (and  M-point  averaging  filter)  satisfies 
(P  (X  <  Y|iT0)  <  §).  With  $  (a-)  as  the  cumulative  distri¬ 
bution  function  of  the  A f  (0, 1)  density  at  the  point  x,  this  re¬ 
sults  in  (V  <  (#)  +  2 yjt  (A0  +  ^L(tissue)))  •  To 

give  a  satisfactory  observation,  we  also  bound  the  miss  prob¬ 
ability  for  detecting  a  ballistic  photon  by  the  same  modified 
value  ( Y)  such  that  (P  (X  >  y\Hi)  <  j^).  Using  the  miss 
bounds,  we  determine  the  lower  bound  on  the  necessary  aver¬ 
aging  window  size  (M)  to  image  a  fixed  K-point  array. 


M  > 


('/(A) 


T  A tissue )  \J (Yo  T  ^tis sue-\-tumor) 


(5) 


The  minimum  width  of  the  occluding  tumor  (wlallistic) 
that  can  be  reliably  resolved  for  a  given  parameterized  turbid 
medium  can  now  be  derived.  Using  the  lower  bound  for  M 
found  in  Eqn  5,  we  can  solve  for  the  lower  bound  on  the  width 

using  ^ballistic  =  y/FOK*M-  Due  to  M  being  a  function 
of  the  tumor  depth  df  ,  we  can  also  numerical  solve  for  the 
minimum  tumor  depth  possible  for  a  parameterized  system. 


2.2.2.  False  Discovery  Rate 

While  the  closed  form  resolution  bounds  for  the  Bonferroni 
Correction  were  derived,  it  is  a  very  conservative  approach 
and  may  obtain  a  poor  reconstruction  in  order  to  avoid  false 
alarm  errors.  To  improve  the  reconstruction,  we  can  increase 
the  resolution  by  decreasing  the  M  value  (averaging  filter 
size),  and  then  use  a  modified  False  Discovery  Rate  (FDR) 
algorithm  [6]  for  the  multiple  point  test.  To  obtain  the  thresh¬ 
old,  take  the  K  number  of  observed  signal  values  and  deter¬ 
mine  the  p-value  (p,)  under  each  observed  value  (X'i). 


Pi 


PixKXWHo)  (6) 

$  (VM  (x\  -  2  *  t  (A0  +  AL(tissiie))^  )  (7) 


To  choose  the  FDR  threshold  (7'),  take  the  threshold  cor¬ 
responding  to  the  largest  index  (n)  such  that  pn  < 

1 

1  —  (1  —  a)A'+1^".  In  practice.  False  Discovery  Rate  will 
result  in  a  less  conservative  reconstruction,  but  it  cannot  be 
analyzed  to  obtain  closed  form  bounds. 

3.  CONVENTIONAL  IMAGING  ANALYSIS 

In  the  conventional  imaging  regime,  there  is  no  time-gating 
mechanism  and  all  the  photons  that  reach  the  detector  over 
a  long  acquisition  time  will  be  observed  (acquisition  time 
^>  (^)  =  direct  line-of-sight  flight  time).  Therefore,  a  large 
number  of  photons  sent  through  the  medium  will  be  collected 
by  the  detector.  A  problem  occurs  here,  too  —  while  the 
signal-to-noise  ratio  is  high  due  to  the  large  number  of  pho¬ 
tons,  the  average  number  of  scattering  events  on  each  pho¬ 
ton  collected  will  also  be  high.  As  the  number  of  scatter¬ 
ing  events  increases  for  a  photon,  the  less  the  photon  will  re¬ 
tain  the  spatial  resolution  of  the  occluding  object.  The  lack 
of  spatial  information  results  in  a  blurred  observation.  Us¬ 
ing  random  walk  theory  in  [2],  it  is  possible  to  solve  for  the 
minimum  width  of  an  occluding  tumor  that  is  reliably  re¬ 
solved  using  the  conventional  imaging  regime.  The  width 
is  found  using  the  photon  mean-time-of-flight  {At},  which 
can  be  numerically  solved  as  a  function  of  the  parameters 
of  the  medium  (/r™,/z™).  The  minimum  width  is  equal  to 


4.  OPTIMAL  RESOLUTION  TRADEOLLS 

Ideally,  one  should  choose  the  imaging  system  (ballistic  or 
conventional)  that  reliably  resolves  the  smallest  possible  ob¬ 
ject  (wt  =  min  wiaUistic)).  The  decision  test  w*om, 

conv 

u’lal ,  using  the  minimum  resolvable  sizes  derived 

ballistic 

above  becomes:  0.408  ( 2  ^  Using 

the  lower  bound  of  M  from  Eqn.  5,  one  can  solve  for  the 
critical  distance  ( dm  =  dcriticai )-  the  maximum  distance  at 
which  ballistic  still  offers  superior  resolution  relative  to  con¬ 
ventional  imaging. 

4.1.  Simulation  Study  -  Breast  Tissue 

We  now  present  a  simulation  study  of  imaging  malignant 
breast  tissue  through  healthy  breast  tissue.  For  this  simu¬ 
lation,  a  10cm  x  10cm  FOV  (as  a  K=2562  point  array)  is 
defined  with  circular  occluding  objects  of  diameter  0.2  cm, 
0.4  cm,  0.8  cm,  2.0  cm,  and  4.0  cm.  From  [3],  we  use  the 
scattering  and  absorption  coefficients  of  the  two  tissue  types. 
Using  a  specified  depth  of  the  occluding  tumor  (d*=0.25cm) 
and  false  alarm  rate  (a  =  0.05),  we  solve  for  the  minimum 
width  of  the  occluding  tumor  under  the  ballistic  observation 


with  Bonferroni  Correction  testing  (wlai(bon)’  al°n8  with  the 
size  of  the  averaging  filter  M),  and  the  conventional  imag¬ 
ing  observations  minimum  occluding  width  (w*onv).  To  ob¬ 
tain  a  better  observation  of  the  model,  a  higher  resolution 

C Wbal(fdr )  =  !U;Lt(6On))0bserVati0nUsingtheFalseDisC0V- 

ery  Rate  (FDR)  method  was  also  simulated.  Three  medium 
distances  ( dm  =  1.6, 1.7, 1.8  cm)  are  used  to  show  the  effects 
that  the  distance  of  the  medium  has  on  the  number  of  ballistic 
photons  received. 


1(A)  Ballistic  Obs.  1(B)  Bonferroni  1(C)  FDR 


•  •  ## 


2(A)  Ballistic  Obs.  2(B)  Bonferroni  2(C)  FDR 

1 


3(A)  Ballistic  Obs.  3(B)  Bonferroni  3(C)  FDR 

Figure  2.  Simulations  under  FOV  =  10x10cm,  d}  =  0.25  cm 
(1)  dr  =  1.6  cm,  (2)  dr  =  1.7  cm,  (3)  dm  =  1.8  cm 


dr 

M 

^ bal(bon) 

^ conv 

pSNRbon 

pSNRfdr 

1.6 

31 

0.217 

1.012 

15.52 

16.34 

1.7 

301 

0.678 

1.04 

11.27 

11.33 

1.8 

3195 

2.21 

1.07 

5.73 

6.21 

Table  1.  Derived  properties  of  the  malignant/healthy 
breast  tissue  environment  (distances  in  cm,  pSNR  in  dB) 


Figure  2  shows  the  effect  of  distance  on  the  resolution 
of  the  observed  image.  One  can  observe  the  ballistic  obser¬ 
vation  with  Bonferroni  Correction  testing  width  size  increas¬ 
ing  dramatically  as  d  — >  dcritiCai  (numerically  found  here  to 
be  =  1.73cm),  this  is  due  to  the  Bonferroni  Correction  being 
a  very  conservative  estimate.  For  d  >  dcriticai  one  can  ob¬ 
serve  that  the  conventional  imaging  observation  performs  sig¬ 
nificantly  better  than  the  ballistic  observation  under  Bonfer¬ 
roni  ( w\al(bon )  >  wtconv )•  As  stated  previous,  the  Bonferroni 
Method  obtains  a  lower  bound  for  the  resolution  of  the  obser¬ 
vation  while  retaining  the  false  alarm  rate.  Using  the  higher 


resolution  False  Discovery  Rate  (FDR)  approach,  we  obtain  a 
higher  SNR  than  the  conservative  Bonferroni  approach  while 
maintaining  the  false  alarm  rate.  A  point  of  interest  is  that  the 
filter  window  will  generally  contain  both  occluding  and  non¬ 
occluding  points  (non-ideal  case),  but  this  does  not  greatly 
degrade  the  reconstruction  quality. 

5.  CONCLUSIONS 

Using  the  decision-theoretic  analysis  approach,  it  is  shown 
that  the  resolution  of  the  turbid  medium,  the  smallest  reli¬ 
ably  resolved  object  (in  terms  of  the  width  w*)  for  both  the 
conventional  imaging  and  ballistic  regimes,  can  be  derived. 
For  the  ballistic  regime,  a  trade-off  between  resolution,  dis¬ 
tance,  laser  intensity  and  confidence  level  was  shown.  The 
ballistic  regime  was  considered  under  two  multiple  hypoth¬ 
esis  test  method,  the  Bonferroni  Correction  and  False  Dis¬ 
covery  Rate.  Under  the  conservative  Bonferroni  Correction, 
the  optimal  choice  between  ballistic  and  conventional  imag¬ 
ing  was  derived  and  can  be  used  to  find  the  best  reconstruction 
technique  for  a  given  system  as  a  function  of  the  parameters 
of  the  medium.  Using  the  False  Discovery  Rate  approach, 
it  was  shown  how  to  obtain  a  higher  resolution  observation 
while  still  maintaining  a  specified  false  alarm  rate. 
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ABSTRACT 

Recently  technological  advances  now  enable  time-gated  ac¬ 
quisitions  of  photons  at  very  fast  rates.  This  can  allow  one  to 
separate  scattered  and  unscattered  photons  in  transillumina¬ 
tion  imaging.  Time-resolved  transillumination  (TRT)  imag¬ 
ing  opens  the  door  to  an  new  type  of  imaging  through  turbid 
(scattering)  media  such  as  soft  tissue  and  fog/smoke,  and  ex¬ 
citing  potential  applications  in  bioimaging  and  surveillance. 
This  paper  proposes  a  novel  Maximum  Likelihood  based  ap¬ 
proach  to  TRT  image  reconstruction. 

Index  Terms —  Image  reconstruction,  Poisson  processes, 
EM  Algorithm,  Time-resolved  Transillumination  Imaging 

1.  TIME-RESOLVED  TRANSILLUMINATION 
IMAGING 

Recently  technological  advances  now  enable  time-gated  ac¬ 
quisitions  of  photons  at  very  fast  rates,  fast  enough  to  sep¬ 
arately  collect  unscattered  (first  arrival)  and  scattered  (later 
arrival)  photons  in  transillumination  imaging  [1],  We  refer 
to  this  technology  as  time-resolved  transillumination  (TRT) 
imaging,  which  is  described  in  more  detail  in  the  next  para¬ 
graph  and  depicted  in  Fig.  1.  TRT  opens  the  door  to  an  new 
type  of  imaging  turbid  (scattering)  media  (e.g.,  soft  tissue, 
fog/smoke).  The  ability  to  separately  detect  unscattered  or 
ballistic  photons  can  enable  much  higher  resolution  imaging 
than  possible  using  conventional  imaging  devices,  and  this 
has  exciting  potential  applications  in  bioimaging  and  surveil¬ 
lance.  However,  the  number  of  ballistic  photons  decays  ex¬ 
ponentially  fast  as  the  thickness/depth  of  the  turbid  medium 
increases.  Therefore,  the  high  resolution  information  that  is 
available  is  also  in  a  very  low  SNR  regime.  This  paper  ex¬ 
plores  a  novel  approach  to  TRT  image  reconstruction. 

In  more  detail,  the  TRT  imaging  problem  involves  of  pho¬ 
tons  traveling  through  a  turbid  medium  from  an  source,  through 
an  object  plane,  and  then  onto  an  imaging  plane  as  depicted 

This  work  was  partially  supported  by  DARPA/AFOSR  grant  number 
FA9550-06- 1-0047. 


in  Fig.  1.  Photons  traveling  through  a  scattering  medium  can 
be  roughly  classified  into  three  types:  ballistic,  snake,  and 
diffuse.  Ballistic  photons  experience  no  scattering  and  travel 
in  a  direct  line  of  sight  arriving  first  at  detectors  in  the  im¬ 
age  plane.  Because  of  the  lack  of  scattering,  ballistic  pho¬ 
tons  retain  their  spatial  information  and  arrive  at  the  imag¬ 
ing  plane  at  the  same  relative  location  as  sent  from  the  ob¬ 
ject  plane.  Snake  photons  experience  some  slight  scattering 
through  the  medium,  this  scattering  will  cause  these  photons 
to  arrive  later  than  the  ballistic  photons  and  likely  in  a  slightly 
different  location  than  sent  from  the  object  plane.  Diffuse 
photons  experience  large  amounts  of  scattering  and  arrive  at 
the  image  plane  having  lost  most  of  their  point  of  origin  infor¬ 
mation.  Due  to  the  large  number  of  scattering  events  through 
the  medium,  the  diffuse  photons  will  travel  the  farthest  dis¬ 
tance  to  the  image  plane,  and  therefore  will  arrive  after  the 
snake  and  ballistic  photons.  While  the  inherent  spatial  in¬ 
formation  decreases  in  order  of  ballistic,  snake  and  diffuse 
photons,  the  number  of  photons  (and  hence  inherent  SNR) 
increases  in  the  same  order.  So  we  are  face  with  high  resolu¬ 
tion,  low  SNR  data  at  one  extreme  (ballstic),  and  low  resolu¬ 
tion,  high  SNR  data  at  the  other  (diffuse).  Furthermore,  the 
diffusion  and  SNR  parameters,  which  characterize  the  under¬ 
lying  point  spread  function  (PSF),  are  not  known  precisely  in 
practice. 


Fig.  1.  Example  of  Photons  Through  a  Scattering  Medium 


The  TRT  image  reconstruction  problem  is  essentially  a 
statistical  inverse  problem  (a  particular  form  of  photon-limited 
image  reconstruction),  but  to  the  best  of  our  knowledge  our 
work  here  is  the  first  to  formally  pose  it  as  such.  Due  to 
uncertainties  in  the  PSF,  the  TRT  problem  somewhat  related 
to  blind  image  deconvolution  (BID).  While  BID  is  a  well- 
studied  problem,  there  are  several  unique  aspects  in  the  TRT 
imaging  problem  that  make  it  quite  different  from  standard 
BID  problems.  In  particular,  the  distinctive  features  of  TRT 
imaging  include  the  photon-limited  nature  of  the  data,  the 
time-gated  data  acquisition  which  in  effect  yields  information 
at  multiple  spatial  resolutions  that  can  be  fairly  well  char¬ 
acterized  via  a  diffusion  equation,  and  most  importantly  the 
availability  of  “unblurred”  data  corresponding  to  the  ballistic 
photons.  For  these  reasons  “off-the-shelf"  BID  algorithms  are 
not  directly  applicable  to  TRT  image  reconstruction. 

This  paper  is  organized  as  follows.  In  Section  2  we  pro¬ 
pose  a  statistical  model  for  TRT  imaging  through  a  homoge¬ 
neous  turbid  medium.  In  Section  3  we  review  a  multiscale 
Poisson  denoising  technique  that  can  be  applied  directly  to 
the  ballistic  data,  and  will  also  be  an  integral  component  of 
our  overall  reconstruction  procedure.  In  Section  4  we  show 
that  an  Expectation-Maximization  (EM)  algorithm  based  on 
the  Poisson  denoising  scheme  can  be  used  to  solve  the  image 
reconstruction  problem  when  one  has  perfect  knowledge  of 
the  scattering  properties  of  the  medium.  In  Section  5  we  pro¬ 
pose  a  novel  scheme  (based  on  the  EM  algorithm)  for  com¬ 
puting  the  joint  Maximum  Penalized  Likelihood  Estimate  of 
the  underlying  image  intensity  and  key  diffusion  and  SNR 
parameters  of  the  scattering  environment.  Section  6  evaluates 
the  performance  of  our  scheme  in  simulations  and  concluding 
remarks  are  made  in  Section  7. 

2.  A  STATISTICAL  MODEL  OF 
TRANSILLUMINATION  IMAGING 

The  basic  statistical  model  we  propose  for  TRT  imaging  through 
homogeneous  turbid  media  is  as  follows.  Assume  that  we 
have  k  time-resovled  “snapshots”  ,  each  acquired  over  dis¬ 
joint  time  intervals  T), . . . ,  Tj~,  with  Y)  denoting  the  ballistic 
time  interval.  Assume  that  these  intervals  form  a  uniform  par¬ 
tition  of  the  overall  observation  interval  T.  Let  X\ , ,  Xk 
denote  the  photon  data  acquired  in  each  interval,  respectively. 
Specifically,  the  data  X,  are  acquired  in  the  form  of  an  n- 
pixel  image,  and  for  our  mathematical  exposition  we  assume 
that  the  columns  of  this  image  are  “stacked”  to  form  an  n  x  1 
vector.  Each  pixel  value  in  A,  is  simply  the  number  of  pho¬ 
tons  detected  at  the  corresponding  location  during  the  time 
interval  I).  Each  image  is  Poisson  distributed  according  to 
the  following  model: 

X,  ~  Poisson(ctjP,;A),  i  =  l,...,k,  (1) 

where  A  denotes  the  underlying  n  x  1  image  intensity  func¬ 
tion,  Pi  denotes  the  n  x  n  photon  transition  matrix  from  the 


emission  (source)  plane  to  the  detection  (image)  plane,  and 
a-i  >  0  is  a  scalar  gain  factor.  The  transition  matrices  are 
functions  of  time  and  a  scalar  diffusion  bandwidth  parameter 
denoted  by  a2.  In  particular,  according  to  the  basic  physics  of 
photon  propagation  through  turbid  media  [1],  row  s  of  Pi  is 
a  probability  mass  function  modeled  by  a  sampled  Gaussian 
density  with  mean  s  and  variance  proportional  to  a2/:,,  where 
ti  denotes  the  midpoint  of  the  ;-th  acquisition  time-interval. 
Thus,  the  transition  matrices  are  parametric  functions  of  the 
form  Pi  =  P(a2ti).  We  assume  that  A  is  normalized  such 
that  the  intensities  in  A  sum  to  one  (i.e.,  YHj= i  At  =  !)•  Also 
the  transition  matrices  are  normalized  so  that  //.,  =  /’A  also 
satisfies  A H,j  =  1-  With  these  normalizations,  it  is  easy 

to  verify  that  the  “total”  intensity  (integrated  spatially  over 
the  image  plane  and  temporally  over  the  i-th  time-interval)  is 
a:,;.  Furthermore,  we  adopt  the  convention  that  the  ballistic 
image  resolution  is  the  finest  spatial  resolution  available  and 
assume  that  Pi  =  InXn ,  the  n  x  n  identity  matrix.  We  also 
assume  that  the  images  acquired  in  each  time-interval  are  sta¬ 
tistically  independent,  and  so  the  joint  distribution  function  of 
the  entire  data  record  X  =  [Xj', . . . ,  X^]T  (the  superscript 
T  denotes  matrix  transposition)  is 

X  ~  Poisson(PA),  (2) 

where  P  is  the  kn  x  n  transition  matrix  obtained  by  stacking 
the  matrices  aqPi, . . . ,  OikPk\  i.e.,  P  =  [ot\Pi  , . . . ,  ctfcPj]T. 

Let  us  constrast  this  imaging  system  with  a  conventional 
system  in  which  the  photon  detections  are  not  time-resolved. 
In  this  case,  we  acquire  the  aggregated  photon  image  Xa  = 
Xi  +  •  •  •  +  Xk  which  obeys  the  model 

Xa  ~  Poisson(PaA),  (3) 

where  Pa  =  2_\=1  aiPi •  We  will  see  that  the  extra  “infor¬ 
mation”  available  in  the  time-resolved  photon  acquisition  can 
significantly  improve  our  ability  to  estimate  the  underlying 
intensity  A. 

3.  PHOTON-LIMITED  IMAGE  DENOISING 

The  ballistic  photon  image  X1  has  high  spatial  resolution  but 
extremely  poor  SNR  due  to  the  very  limited  number  of  bal¬ 
listic  photons.  As  a  starting  point  for  our  work,  let  us  con¬ 
sider  the  problem  of  estimating  the  underlying  image  inten¬ 
sity  based  on  the  ballistic  photon  data  alone.  This  boils  down 
to  a  Poisson  image  “denoising”  problem,  which  has  recently 
received  a  significant  amount  of  attention  in  the  image  pro¬ 
cessing  and  statistics  literature. 

One  state-of-the-art  Poisson  denoising  scheme  is  based 
on  the  recursive  dyadic  partition  (RDP)  framework  proposed 
in  [2],  This  scheme  is  a  Poisson  analog  of  the  more  fa¬ 
miliar  wavelet  denoising  methods  developed  for  the  classi¬ 
cal  “signal+noise”  model.  Also  like  wavelet  denoising,  addi¬ 
tional  improvements  in  denoising  quality  are  possible  using 


a  translation-invariant  version  of  the  basic  approach  [3,  4], 
which  can  be  computed  in  0(n  log  n)  operations. 

The  Poisson  denoising  method  will  be  an  integral  compo¬ 
nent  of  our  EM  algorithm  for  time-resolved  transillumination 
imaging.  The  EM  algorithm  optimally  combines  information 
from  the  entire  record  time-resolved  photons  (i.e.,  from  the 
ballistic,  quasidiffuse  and  diffuse  limits).  But  before  moving 
on,  let  us  illustrate  the  performance  of  the  method  by  esti¬ 
mating  the  underlying  image  intensity  using  only  the  ballistic 
photon  data.  Figure  2  depicts  the  results  of  denoising  a  bal¬ 
listic  photon  image  using  the  methodology  described  above 
(specifically,  we  employ  the  translation-invariant  Haar  esti¬ 
mation  scheme  described  in  [4]  and  implemented  in  the  su¬ 
perb  Matlab  package  developed  by  Prof.  Rebecca  Willett  [5]). 
Note  that  the  denoising  method  results  in  an  intensity  estimate 
that  reduces  “noise”  while  preserving  edges  and  other  details. 


5.  ADAPTING  TO  UNKNOWN  TURBID  MEDIA 

One  of  the  major  challenges  in  practical  imaging  problems  is 
that  the  characteristics  of  the  turbid  medium  are  usually  not 
known  precisely.  In  particular,  the  values  of  system  param¬ 
eters  {cti}  and  a2  are  unknown  and  therefore  must  be  esti¬ 
mated  along  with  A.  We  can  formulate  this  as  a  joint  MLE 
problem,  seeking  to  find  values  of  the  system  parameters  and 
A  which  jointly  maximize  the  Poisson  likelihood  function  (or 
penalized  likelihood).  At  first  glance,  it  may  appear  that  this 
joint  MLE  problem  may  be  intractable,  but  it  turns  out  to  have 
a  rather  simple  solution  which  is  one  of  the  main  contributions 
of  this  paper. 

The  basic  solution  approach  is  as  follows.  First,  the  MLEs 
of  the  gain  factors  {ct,}i=i  can  be  computed  separately  from 
A  and  a2  due  to  the  following  observation.  Consider  the 
statistics 

n 

Si  =  ^ Xiyj ,  i  =  l,...,k, 
j= i 

the  subscript  j  indexes  the  pixels  in  each  image.  These  statis¬ 
tics  are  simply  the  total  photon  counts  for  each  image.  Due  to 
the  normalization  of  A  and  the  matrices  { /  ' }  in  our  model,  it 
follows  that 


Fig.  2.  Example  of  ballistic  photon  image  MPLE  denoising. 


4.  AN  EM  ALGORITHM  FOR  IMAGE 
RECONSTRUCTION 

If  the  attenuation  factors  {c^}  and  scattering  matrices  {Pi} 
are  known,  then  the  MLE  of  A  given  X  (or  given  Xa)  can  be 
computed  using  the  well-known  Expectation-Maximization 
(EM)  algorithm  [6].  For  the  general  problem,  the  standard 
E  and  M  steps  are  computed  using  X  and  P,  while  in  the  case 
of  aggregated  photon  data  the  steps  employ  Xa  and  Pa.  The 
algorithm  is  initialized  with  a  starting  guess  for  MLE  of  A 
(e.g.,  all  pixels  set  to  1). 

Under  mild  conditions,  the  sequence  of  estimates  con¬ 
verges  to  an  MLE  solution.  Unfortunately,  because  P  is  usu¬ 
ally  poorly-conditioned,  the  exact  MLE  solution  is  usually 
undesirable.  For  example,  in  the  case  of  the  ballistic  data 
alone,  the  MLE  of  A  is  simply  X\,  which  as  seen  from  Fig¬ 
ure  2,  is  highly  variable  and  typically  has  a  very  poor  MSE. 
So,  instead  of  seeking  the  MLE  we  aim  to  recover  a  Maxi¬ 
mum  Penalized  Likelihood  Estimate  (MPLE),  using  the  Pois¬ 
son  denoising  criterion  described  in  Section  3  as  our  penalty 
term.  This  MPLE  approach  was  first  proposed  in  [7].  The 
MPLE  can  also  be  computed  using  the  EM  algorithm.  In  this 
case,  the  E-Step  remains  the  same  as  above  and  the  M-Step  is 
computed  by  applying  the  translation-invariant  denoising  al¬ 
gorithm  [4,  5]  to  the  usual  M-step  result  prior  to  re -computing 
the  E-Step  (see  [7]  for  further  details). 


Si  ~  Poisson(ct;),  i  =  1, . . .  ,k. 

It  is  well-known  that  the  conditional  distribution  of  X,  given 
Si  is  multinomial  with  parameters  //,  =  PjA  (see  [2]).  There¬ 
fore,  the  likelihood  factorizes  into  Poisson  factors,  each  in¬ 
volving  one  pair  (cti,  Si),  and  multinomial  factors,  each  in¬ 
volving  A  and  one  triple  (Pi,  Xi,  Si).  Consequently,  the 
MLEs  of  the  gain  factors  {ct;}  can  be  obtained  separately 
from  A  and  a2  and  are  given  by  the  simple  formula  a,  =  ,S',, 
i  =  1 ,...  ,k.  Now  recall  that  the  matrices  {P,}  are  para¬ 
metric  in  a2.  To  find  the  MLEs  of  A  and  a2  we  consider 
a  range  of  candidate  values  for  a 2  and  for  each  one  we  use 
the  EM  algorithm  described  above  (with  each  a,  set  to  its 
MLE  ai)  to  compute  the  MPLE,  denoted  by  A(cr2),  and  the 
corresponding  maximum  penalized  likelihood  value,  denoted 
L(a2).  This  can  be  done  exhaustively  (over  a  discretized 
range  of  a2  values)  or  systematically  (assuming  unimodality 
of  the  maximum  penalized  likelihood  as  a  function  of  a2  and 
employing  a  bisection  method).  The  value  of  a2  that  results 
in  the  highest  penalized  likelihood  value  L(cr2)  is  the  MLE  of 
a2,  denoted  by  a2.  The  MLE  of  A  is  then  A(cr2). 

6.  A  TRT  IMAGING  EXAMPLE 

The  potential  of  the  proposed  EM  algorithm  is  evaluated  in 
the  following  simulation  study.  Using  the  A  intensity  function 
depicted  in  Figure  2a,  a  ballistic  image  is  generated  using  low 
photon  count  Poisson  data  and  a  diffuse  image  is  generated  by 
first  blurring  A  using  a  high  variance  blurring  kernel  and  then 


generating  a  large  number  of  photons  from  the  blurred  inten¬ 
sity  function.  The  reconstruction  based  only  on  the  ballistic 
photon  data  is  depicted  in  Figure  2c,  and  the  reconstruction 
using  only  the  diffuse  image  and  assuming  perfect  knowl¬ 
edge  of  the  blurring  variance  cr2  is  depicted  in  Figure  3a). 
Figure  3b  shows  the  reconstruction  based  on  both  ballistic 
and  diffuse  photon  data  and  assuming  perfect  knowledge  of 
attenuation  factors  (ai,  0:2)  and  blurring  variance  (cr),  using 
the  EM  algorithm  described  in  Section  4.  Using  the  adap¬ 
tive  MLE  scheme  described  in  Section  5,  the  attuenation  fac¬ 
tors  and  blurring  variance  were  jointly  estimated  along  with 
A,  and  the  resulting  intensity  estimate  shown  in  (Figure  3c). 
Table  1  summarizes  the  reconstuction  errors  over  10  indepen¬ 
dent  trials  of  the  experiment.  As  expected,  the  combined  (bal- 
listic+diffuse)  oracle  reconstruction  performs  best,  with  the 
combined  reconstruction  with  ML  estimation  of  the  gain  and 
diffusion  parameters  close  behind. 

The  algorithm  developed  in  Section  5  the  final  reconstruc¬ 
tion  is  the  estimate  with  largest  likelihood  value  chosen  from 
a  set  of  reconstructed  images  generated  using  different  cr2 
values.  In  all  our  experiments,  the  likelihood  appears  to  be 
unimodal  in  cr2  (one  such  likelihood  is  seen  in  (Figure  3d)), 
allowing  for  systematic  searches  such  as  a  bisection  method. 
We  conjecture  that  the  likelihood  is  always  unimodal  in  cr2, 
but  we  have  not  yet  proved  this  result. 


(a)  MPLE  Diffuse 


(b)  MPLE  Ballistic  +  Diffuse 


(d)  Log -Likelihood  Plot 

Fig.  3.  Examples  of  TRT  image  reconstruction. 


(c)  MPLE  Ballistic  +  Diffuse 
w/  Unknown  parameters 


Table  1.  Summary  of  TRT  image  reconstruction  errors. 


Image 

PSNR  (dB) 

Ballistic  Image 

4.08 

Diffuse  Image 

7.00 

MPLE  Ballistic 

11.90 

MPLE  Diffuse 

9.42 

MPLE  Ballistic  +  Diffuse 

12.92 

MPLE  Ballistic  +  Diffuse  w /  unknown  params 

12.71 

7.  CONCLUSIONS 

This  paper  proposed  a  novel  MLE  reconstruction  algorithm 
for  TMT  imaging.  The  algorithm  is  based  on  a  combination 
of  the  EM  algorithm  and  Poisson  denoising.  Our  simulation 
study  demonstrates  the  potential  of  our  approach,  in  particular 
indicating  the  added  benefit  of  optimally  fusing  ballistic  and 
diffuse  photon  data.  Our  future  work  includes  the  application 
of  this  theory  to  real-data  experiments,  detailed  analysis  of 
fundamental  performance  limits  in  TRT  imaging,  and  exten¬ 
sions  to  inhomogeneous  media. 
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